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1 KINEMATICS OF MOVING FRAMES 

1.1 Rotation of Reference Frames 

We denote through a subscript the specific reference system of a vector. Let a vector ex- 
pressed in the inertial frame be denoted as x, and in a body-reference frame Xb- For the 
moment, we assume that the origins of these frames are coincident, but that the body 
frame has a different angular orientation. The angular orientation has several well-known 
descriptions, including the Euler angles and the Euler parameters (quaternions). The former 
method involves successive rotations about the principle axes, and has a solid link with the 
intuitive notions of roll, pitch, and yaw. One of the problems with Euler angles is that for 
certain specific values the transformation exhibits discontinuities. Quaternions present a 
more elegant and robust method, but with more abstraction. We will develop the equations 
of motion using Euler angles. 

Tape three pencils together to form a right-handed three-dimensional coordinate system. 
Successively rotating the system about three of its own principal axes, it is easy to see 
that any possible orientation can be achieved. For example, consider the sequence of [yaw, 
pitch, roll]: starting from an orientation identical to some inertial frame, rotate the movable 
system about its yaw axis, then about the new pitch axis, then about the newer still roll 
axis. Needless to say, there are many valid Euler angle rotation sets possible to reach a given 
orientation; some of them might use the same axis twice. 




Figure 1: Successive application of three Euler angles transforms the original coordinate 
frame into an arbitrary orientation. 



A first question is: what is the coordinate of a point fixed in inertial space, referenced to 
a rotated body frame? The transformation takes the form of a 3x3 matrix, which we now 
derive through successive rotations of the three Euler angles. Before the first rotation, the 
body-referenced coordinate matches that of the inertial frame: Now rotate the 

movable frame yaw axis (z) through an angle </>. We have 



Xu 



COS ( 

- sin i 



sin (f) 
cos (j) 
1 



R(<f>)x 



6 • 



(1) 
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1 KINEMATICS OF MOVING FRAMES 



Rotation about the z-axis does not change the ^-coordinate of the point; the other axes are 
modified according to basic trigonometry. Now apply the second rotation, pitch about the 
new y-axis by the angle 9: 



^2 



X 



cos 9 
1 
sin^ 



- sin 9 


cos 9 



Xl 



R(9)x 



b ■ 



(2) 



Finally, rotate the body system an angle ip about its newest x-axis: 



^3 



x 



1 







cos ip sin tp 
- sin ip cos ip 



x, 



R(tP)x 



^2 
b ■ 



(3) 



This represents the location of the original point, in the fully-transformed body-reference 
frame, i.e., x^. We will use the notation Xf, instead of x^ from here on. The three independent 
rotations can be cascaded through matrix multiplication (order matters!): 



x b = R{i;)R{9)R((P)x 

c9ccf) 

—cipstj) + sif)s9c<p 
sipscj) + af)s9c(fi 

R{(p,9,i;)x. 



(4) 



c9s(f) —s9 
ci/jc<j) + sips9scj) siJjc9 
-sipc(f) + cips9s(f> apc9 



X 



All of the transformation matrices, including R(<p,9,ip), are orthonormal: their inverse is 
equivalent to their transpose. Additionally, we should note that the rotation matrix R 
is universal to all representations of orientation, including quaternions. The roles of the 
trigonometric functions, as written, are specific to Euler angles, and to the order in which 
we performed the rotations. 

In the case that the movable (body) reference frame has a different origin than the inertial 
frame, we have 



x = x + R T Xb, 

where x is the location of the moving origin, expressed in inertial coordinates. 



(5) 



1.2 Differential Rotations 

Now consider small rotations from one frame to another; using the small angle assumption 
to ignore higher-order terms gives 



R ~ 



59 



1 

-Sip 



-59 

5i) 
1 



(6) 



1.2 Differential Rotations 



3 



Sep -89 
■6<f> 5i> 
S9 -Sip 



3x3) 



where ^3x3 donotes the identity matrix. R comprises the identity plus a part equal to the 
(negative) cross-product operator [— SEx], where SE = [Sip, S9, Sep], the vector of Euler 
angles ordered with the axes [x,y, z]. Small rotations are completely decoupled; the order of 
the small rotations does not matter. Since R^ 1 = R T , we have also R^ 1 = / 3X 3 + SEx; 



x b = x — SE x x (7) 
x = Xb + SE x x b . (8) 

We now fix the point of interest on the body, instead of in inertial space, calling its location 
in the body frame r (radius). The differential rotations occur over a time step St, so that 
we can write the location of the point before and after the rotation, with respect to the first 
frame as follows: 



x(t) = f (9) 
x(t + St) = R T r = f+SEx f. 

Dividing by the differential time step gives 

Sx 
It 

where the rotation rate vector uj ~ dE/dt because the Euler angles for this infinitesimal 
rotation are small and decoupled. This same cross-product relationship can be derived in 
the second frame as well: 



SE 



x r 



St 

uj x f, 



(10) 



x b {t) 
x b (t + St) 



Rf 



r. 



SExf 



:n) 



such that 



Sx b SE _ 

w = j x j (12) 

= UJ X f, 

On a rotating body whose origin point is fixed, the time rate of change of a constant radius 
vector is the cross-product of the rotation rate vector uj and the radius vector itself. The 
resultant derivative is in the moving body frame. 
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1 KINEMATICS OF MOVING FRAMES 



In the case that the radius vector changes with respect to the body frame, we need an 
additional term: 



dx b 
~dT 



dr 

uxr + —. 
at 



Finally, allowing the origin to move as well gives 



dxb 
~dt 



dr dx Q 



(13) 



(14) 



This result is often written in terms of body-referenced velocity v: 

dr 

v = u}xf+ — + v , (15) 

where v is the body-referenced velocity of the origin. The total velocity of the particle is 
equal to the velocity of the reference frame origin, plus a component due to rotation of this 
frame. The velocity equation can be generalized to any body-referenced vector /: 



dj_ 
dt 



dt 



+ x /. 



(16) 



1.3 Rate of Change of Euler Angles 

Only for the case of infinitesimal Euler angles is it true that the time rate of change of 
the Euler angles equals the body-referenced rotation rate. For example, with the sequence 
[yaw,pitch,roll], the Euler yaw angle (applied first) is definitely not about the final body yaw 
axis; the pitch and roll rotations moved the axis. An important part of any simulation is 
the evolution of the Euler angles. Since the physics determine rotation rate u, we seek a 
mapping u — > dE/dt. 

The idea is to consider small changes in each Euler angle, and determine the effects on the 
rotation vector. The first Euler angle undergoes two additional rotations, the second angle 
one rotation, and the final Euler angle no additional rotations: 



u = R{ip)R{9) 



+ R(ip) < 



1 -sin^ 

cos i(j sin ip cos 9 
— sin -0 cos ip cos 9 





1 




< 


de 


> 




dt 








> 






( dip 


> 




) % 


















\ dt 


> 



Taking the inverse gives 



dip 

dt 





(17) 



dE 
~dt 



1 sin -0 tan 9 cos ip tan 6 

cos -0 — sin ^ 

sin-0/cos^ cos-0/cos6> 

= T(E)u. 



(18) 



1.4 Dead Reckoning 
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Singularities exist in T at 6 = {n/2, 37r/2}, because of the division by cos#, and hence this 
otherwise useful equation for propagating the angular orientation of a body fails when the 
vehicle rotates about the intermediate y-axis by ninety degrees. In applications where this 
is a real possibility, for example in orbiting satellites and robotic arms, quaternions provide 
a seamless mapping. For most ocean vessels, the singularity is acceptable, as long as it is 
not on the yaw axis! 



1.4 Dead Reckoning 

The measurement of heading and longitudinal speed gives rise to one of the oldest methods 
of navigation: dead reckoning. Quite simply, if the estimated longitudinal speed over ground 
is U, and the estimated heading is (f>, ignoring the lateral velocity leads to the evolution of 
Cartesian coordinates: 



x = U cos (f) 
y = U sin 0. 

Needless to say, currents and vehicle sideslip will cause this to be in error. Nonetheless, some 
of the most remarkable feats of navigation in history have depended on dead reckoning. 



2 VESSEL INERTIAL DYNAMICS 

We consider the rigid body dynamics with a coordinate system affixed on the body. A 
common frame for ships, submarines, and other marine vehicles has the body-referenced x- 
axis forward, y-axis to port (left), and z-axis up. This will be the sense of our body-referenced 
coordinate system here. 



2.1 Momentum of a Particle 

Since the body moves with respect to an inertial frame, dynamics expressed in the body- 
referenced frame need extra attention. First, linear momentum for a particle obeys the 
equality 

F = -j- (mv) (19) 

A rigid body consists of a large number of these small particles, which can be indexed. The 
summations we use below can be generalized to integrals quite easily. We have 

Pl + ^ l = Jt ' (20) 

— * — # 

where Ft is the external force acting on the particle and E4 is the net force exerted by all the 
other surrounding particles (internal forces). Since the collection of particles is not driven 
apart by the internal forces, we must have equal and opposite internal forces such that 
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2 VESSEL INERTIAL DYNAMICS 



N 



££ = o. 



(21) 



Then summing up all the particle momentum equations gives 



N - N d 



(22) 



i=l i=l 



Note that the particle velocities are not independent, because the particles are rigidly at- 
tached. 

Now consider a body reference frame, with origin 0, in which the particle i resides at body- 
referenced radius vector r; the body translates and rotates, and we now consider how the 
momentum equation depends on this motion. 



Figure 2: Convention for the body-referenced coordinate system on a vessel: x is forward, 
y is sway to the left, and z is heave upwards. Looking forward from the vessel bridge, roll 
about the x axis is positive counterclockwise, pitch about the y-axis is positive bow-down, 
and yaw about the z-axis is positive turning left. 

2.2 Linear Momentum in a Moving Frame 

The expression for total velocity may be inserted into the summed linear momentum equation 
to give 



where m = J2i=i m i, an d Vi = v a + u x r im Further defining the center of gravity vector fa 
such that 





(23) 




rriiri 



2.3 Example: Mass on a String 
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N 



i=l 



we have 



N 



3 d v» , d i - — * \ 



(24) 



(25) 



Using the expansion for total derivative again, the complete vector equation in body coor- 
dinates is 



F = 2_^N = m ^— + u x u c + — x r G + to x (u x r G ) 



(26) 



Now we list some conventions that will be used from here on: 



v = {u, v, w} (body-referenced velocity) 

r G = {x G , Vgi z g} (body-referenced location of center of mass) 
^ = {p, 9; r } (rotation vector, in body coordinates) 
F = {X, Y, Z} (external force, body coordinates) . 
The last term in the previous equation simplifies using the vector triple product identity 

uj x (uj x fa) = (uj ■ r G )u) — (u ■ io)f G , 
and the resulting three linear momentum equations are 



X = m 
Y = m 
Z = m 



— + qw-rv+ —z G - —y G + (qy G + rz G )p - (q 2 + r 2 )x G 



(27) 



dv 



dr 



+ ru-pw + —x G 
at dt 



dp 
~dt 



z G + (rz G + px G )q - (r 2 + p 2 )y G 



dw dp dq . . 2 2 . 

— + pv - qu + — y G - —x G + (px G + qy G )r - (p + q )z G 



Note that about half of the terms here are due to the mass center being in a different location 
than the reference frame origin, i.e., f G ^ 0. 

2.3 Example: Mass on a String 

Consider a mass on a string, being swung around around in a circle at speed U, with radius r. 
The centrifugal force can be computed in at least three different ways. The vector equation 
at the start is 



^ / dv Q _ _ duj _ _ ,_ _ . \ 
F = m\ — +ujxv + — xr G + wx(wxr G ) . 
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2 VESSEL INERTIAL DYNAMICS 



2.3.1 Moving Frame Affixed to Mass 

Affixing a reference frame on the mass, with the local x oriented forward and y inward 
towards the circle center, gives 

v = {U,0,0} T 
u = {0,0, U/r} T 
r G = {0,0, 0} T 

such that 

F = muj x v Q = m{0, U 2 /r, 0} T . 
The force of the string pulls in on the mass to create the circular motion. 

2.3.2 Rotating Frame Attached to Pivot Point 

Affixing the moving reference frame to the pivot point of the string, with the same orientation 
as above but allowing it to rotate with the string, we have 

v = {0,0, 0} T 
u = {0,0, U/r} T 
f G = {0,r,0} T 

{0,0, 0}^ 



dt 

dou 

~dt 



{0,0, 0} T , 



giving the same result: 



F = muj x (d? x f G ) = m{0, U 2 /r, 0} rj 



2.3.3 Stationary Frame 

A frame fixed in inertial space, and momentarily coincident with the frame on the mass 
(2.3.1), can also be used for the calculation. In this case, as the string travels through a 
small arc Sip, vector subtraction gives 



Sv = {0,f/sin#,0} T ~ {0,U5ip,0} T . 
Since ip — U/r, it follows easily that in the fixed frame dv/dt = {0, U 2 /r, 0} T , as before. 



2.4 Angular Momentum 
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2.4 Angular Momentum 

For angular momentum, the summed particle equation is 



N - N d 

^2(Mi + f i xF i )=^2r i x —(rriiVi), 

i=i i=i dt 



(28) 



where Mj is an external moment on the particle i. Similar to the case for linear momentum, 
summed internal moments cancel. We have 



N i N 



dv 
dt 



+ UJ X v 



N 



' dtJ 



J2 m i?i x (w x (w X fj)). 

i=l 



The summation in the first term of the right-hand side is recognized simply as mr G , and the 
first term becomes 



mro x 

The second term expands as (using the triple product) 



— + ujxv 



(29) 



N 

x 

1=1 



' du 

~dt 



£*Li m» ((y 2 + ^ 2 )p - (j/tf + Zif)xi) 
= ' E 4 =i m {{xf + - (Zip + Zif)yi) 
k Eili ((a; 2 + % 2 )^ - W + 



Employing the definitions of moments of inertia, 



(30) 



I = 



<-xy 





Ixy 




lyx 




lyz 








Izx 


Izy 


Izz 



N 

t=l 
N 

i=l 
N 

i=l 

N 

lyx ^ ] TTiiXiVi 

i=l 



(inertia matrix) 



(cross-inertia) 
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2 VESSEL INERTIAL DYNAMICS 



N 

i=\ 
N 



<-yz 



i=l 



the second term of the angular momentum right-hand side collapses neatly into Iduj/dt. The 
third term can be worked out along the same lines, but offers no similar condensation: 



N 



N 



mn x ((uj • fi)u - (u ■ u)fi) = rriifi x u(uj ■ fj) 



(31) 



t=i 



£*Li mi(yir - Ziq)(xip + y& + zp) 
Eili mi(zip - Xir)(xip + y t q + ^r) 
Y,f=imi(x i q-yip)(x i p + y i q + Zir) j 

I yz {q 2 - r 2 ) + 7 X2 pg - I xy pr 
I X z(r 2 - p 2 ) + I xy rq - I yz pq } + 
I X y{p 2 ~q 2 ) + lyzPr - I xz qr 




Letting M = {K,M,N} be the total moment acting on the body, i.e., the left side of 
Equation 28, the complete moment equations are 



K 



1-xxP IxyQ Ixz'i' 

{hz - I yy )rq + I yz (q 2 - r 2 ) + I xz pq - I xy pr + 
m [yc(w + pv — qu) — z G (v + ru — pw)} 



(32) 



M 



r yz qp + 



N 



yyq ' v z ' 
. , , I zz )pr + I xz (r 2 -p 2 ) + I xy qr - I y 

m [zg(u + qw — rv) — x G {w + pv — qu)] 



hxP + I zy q + hzf + 

{Iyy - I XX )Pq + I X y(p 2 ~ q 2 ) + IyzPT ~ IxzQT + 

m [xq(v + ru — pw) — y G (u + qw — rv)] . 



2.5 Example: Spinning Book 

Consider a homogeneous rectangular block with I xx < I yy < I zz and all off-diagonal moments 
of inertia are zero. The linearized angular momentum equations, with no external forces or 
moments, are 



2.5 Example: Spinning Book 
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Xx dt + Kzz 


- hvYq 


= 


I ^- + (1 

dt 


- Izz)vr 


= 


dr 

zz di + ' 


- i xx )qp ~- 


= 0. 



We consider in turn the stability of rotations about each of the main axes, with constant 
angular rate fi. The interesting result is that rotations about the x and z axes are stable, 
while rotation about the y axis is not. This is easily demonstrated experimentally with a 
book or a tennis racket. 



2.5.1 x-axis 

In the case of the x-axis, p = Q + dp, q = Sq, and r = <5r, where the 5 prefix indicates a small 
value compared to fl. The first equation above is uncoupled from the others, and indicates 
no change in Sp, since the small term SqSr can be ignored. Differentiate the second equation 
to obtain 



d 2 Sq 9fr 

yy~0fT ^ ' xx ~ ' ~Qt = 

Substitution of this result into the third equation yields 



d 2 Sq 2 

lyylzz ^ {.^xx ^-zz)ij-xx lyy)^ 1 &q 0. 

A simpler expression is 5q + ot5q = 0, which has response Sq(t) = #g(0)e v/3 "*, when 8q(0) = 0. 
For spin about the x-axis, both coefficients of the differential equation are positive, and 
hence a > 0. The imaginary exponent indicates that the solution is of the form Sq(t) = 
5q(0)cosy/at, that is, it oscillates but does not grow. Since the perturbation Sr is coupled, 
it too oscillates. 



2.5.2 y-axis 

Now suppose q = fl+Sq: differentiate the first equation and substitute into the third equation 
to obtain 



d 2 Sp 2 

Izzlxx ^ (lyy ^-xx){Iyy Izz)^ $P 0. 

Here the second coefficient has negative sign, and therefore a < 0. The exponent is real now, 
and the solution grows without bound, following 8p(t) = £p(0)e v/ ~"*. 
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2.5.3 z-axis 

Finally, let r = Q+Sr: differentiate the first equation and substitute into the second equation 
to obtain 



d 2 Sp 2 

lyylxx ^ \Ixx Izz)\Iyy Izz)Q $P 0. 

The coefficients are positive, so bounded oscillations occur. 
2.6 Parallel Axis Theorem 

Often, the mass center of an body is at a different location than a more convenient measure- 
ment point, the geometric center of a vessel for example. The parallel axis theorem allows 
one to translate the mass moments of inertia referenced to the mass center into another 
frame with parallel orientation, and vice versa. Sometimes a translation of coordinates to 
the mass center will make the cross-inertial terms I xy , I yz , I xz small enough that they can be 
ignored; in this case r G = also, so that the equations of motion are significantly reduced, 
as in the spinning book example. 
The formulas are: 



Ixx 


I XX 


+ m(Sy 2 + 5z 2 ) 


lyy 


= lyy 


+ m(5x 2 + 5z 2 ) 


J-zz 


Izz 


+ m(5x 2 + Sy 2 ) 


Iyz 


— Iyz 


— mSydz 


Ixz 


Ixz 


— m5x5z 


Ixy 


— Ixy 


— mSxSy, 



where / represents an MMOI in the axes of the mass center, and 5x, for example, is the 
translation of the x-axis to the new frame. Note that translation of MMOI using the parallel 
axis theorem must be either to or from a frame resting exactly at the center of gravity. 



2.7 Basis for Simulation 

Except for external forces and moments F and M, we now have the necessary terms for 
writing a full nonlinear simulation of a rigid body, in body coordinates. There are twelve 
states, comprising the following components: 

• v , the vector of body-referenced velocities. 

• lo, body rotation rate vector. 

• x, location of the body origin, in inertial space. 
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• E, Euler angle vector. 

The derivatives of body-referenced velocity and rotation rate come from Equations 27 and 32, 
with some coupling which generally requires a 6 x 6 matrix inverse. The Cartesian position 
propagates according to 

'x = R T (E)v 0l (34) 

while the Euler angles follow: 

E = T(E)u. (35) 



3 NONLINEAR COEFFICIENTS IN DETAIL 

The method of hydrodynamic coefficients is a somewhat blind series expansion of the fluid 
force in an attempt to provide a framework on which to base experiments and calculations 
to evaluate these terms. The basic dificulty, i.e. the intractability of the governing equations 
of motion of viscous fluid prohibits, at least today and in the near future, a computation of 
these forces. 

Still, we are not totally ignorant about these forces, since a number of symmetries and basic 
laws can be applied to reduce the number of unknown coefficients. This is the purpose of 
this section. 

The basic assumption in using the method of hydrodynamic coefficients is that the forces 
have no memory effects, i.e. past motions have no impact on the fluid forces at the present 
moment. This is not correct when the flow separates, or when large vortices are shed, 
because then the vorticity in the fluid affects the fluid forces for a considerable time after 
they have been shed - i.e. until they move sufficiently far away from the body. We will show 
later some methods which allow us to incorporate the effect of shed vorticity, because under 
certain conditions such effects can not be ignored. 

We employ the following basic facts and assumptions to derive the fluid forces acting on a 
ship, submarine or vehicle: 

1. We retain only first order acceleration terms. Based on Newton's second law, we expect 
the inertia terms from the fluid to be linearly dependent on acceleration. 

2. We do not include terms coupling velocities and accelerations. Again, based on New- 
ton's second law, we expect inertia forces to depend on acceleration alone. 

3. We consider port /starboard symmetry. Unless there is a reason not to, this is a useful 
property to us in order to eliminate a certain number of coefficients which are either 
zero or very small. In ships, a propeller introduces an asymmetry port /starboard since 
it rotates in a certain direction (unless it is a twin-screw ship with counter-rotating 
propellers), but in such cases we limit this asymmetry to only propeller- related terms. 
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4. We retain up to third order terms. This is a practical consideration and has been 
found to, in general, serve well the purpose of deriving equations which are sufficiently 
accurate in a wide parametric range. It does not constitute an absolute rule, but most 
existing models employ this assumption 

Finally, we find sometimes necessary to use coefficients such as Y\ v \ v , providing a drag-related 
term, whose strict definition would be: 



Y \v\v = T^l( v = 0) (36) 



2 Y 
dvd\v\ 

A Taylor series expansion would not include such terms, so their inclusion is motivated by 
physical arguments. 

3.1 Helpful Facts 

To exploit symmetries, we consider the following simple facts: 

• A function f(x) which is symmetric in x has zero odd derivatives with respect to x at 
x = 0. The proof is to consider the Taylor series expansion of f(x) about x = 0: 

/W = /(0 ) + |(0), + ig(0y + ig(0y + ... (37) 
Symmetry gives that for all x: 



f(x) = f(-x) (38) 

The only way that (38) is satisfied given the expression (37) is that all odd derivatives 
of f(x) must be zero, i.e., for n odd: 

^(O) = (39) 

• A function f(x) which is anti-symmetric in x has zero even derivatives with respect to 
x at x = 0. The proof is exactly analogous to the result above, using the anti-symmetry 
condition: 



f(x) = -f(-x) (40) 

to find that, for n even: 



(41) 



3.2 Nonlinear Equations in the Horizontal Plane 
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3.2 Nonlinear Equations in the Horizontal Plane 

To demonstrate the methodology, we will derive the governing nonlinear equations of motion 
in the horizontal plane (surge, sway and yaw), employing the assumptions above. We will 
not include drag related terms of the form of equation 36 for the time being. We start by 
re-stating the inertia terms driven by external forces, which include the fluid forces: 



m(u — vr — x G r 2 ) = X 
m(v + ur + xcr) = Y 
I zz r + mx G (v + ur) = N (42) 

where we have assumed that w = 0, p = 0, q = 0, yc = 0, z G = 0. 
3.2.1 Fluid Force X 

By denoting the rudder angle as 5, we derive the following expression for the fluid force X, 
valid up to third order: 



X = X e + Xuti + X u u + X uu u 2 + X uuu u 3 + X vv v 2 + X rr r 2 + X SS S 2 + 
X rv rv + X rS r5 + X vS vS + X vvu v 2 u + +X rru r 2 u + X S 6 u S 2 u + 
X r s u rSu + X rvu rvu + +X v $ u v8u + X rv $rv8 (43) 

We have used three basic properties, i.e., that the fluid force X, independent of the forward 
velocity, must be: 

1. a symmetric function of v when r = and 5 = 0; 

2. a symmetric function of r when v = and 5 = 0; 

3. a symmetric function of 5 when r = and v = 0; 

This is a result of port /starboard symmetry and is expressed as: 



X(u, v, r = 0, 5 = 0) = X(u, -v,r = 0,5 = 0) (44) 
X(u, v = 0,r,5 = 0)= X(u, v = 0,-r,S = 0) (45) 
X(u, v = 0, r = 0, 5) = X(u, v = 0,r = 0, -5) (46) 

These relations imply that all odd derivatives of X with respect to v at v = are zero, when 
r = and 5 = 0; and similarly for r and 5. For example: 



(u,v = 0,r = 0,5 = 0) = 



(47) 
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which implies that X vvv = 0. Also, since relations such as (44) hold for all forward velocities 
u, it means that derivatives with respect to u of expression (44) are aso true. For example, 
from (47) we derive: 



d 4 X 

' ' (u,v = 0,r = 0,5 = 0) = (48) 



dv 3 du 



or equivalently X vvvu = 0. 

In summary, the symmetries provide the following zero coefficients: 



X v 


= 0; 




= 0; 


X vu 


= 0; 


X vuu 


X r 


= 0; 


Y 


= 0; 


Y 

^-ru 


= 0; 


X ruu 


X s 


= 0; 


Xsss 


= 0; 


Xsu 


= 0; 


Xsuu = 



3.2.2 Fluid Force Y 

In the case of the fluid foce Y, the symmetry implies that this force must be an antisymmetric 
function of v when r = 0, S = 0; and likewise for r and S, i.e.: 



Y(u,v,r = 0,5 = 0) = -Y(u,-v,r = 0,6 = 0) 
Y(u,v = 0,r,5 = 0) = -Y(u,v = 0,-r,5 = 0) 
Y(u,v = 0,r = 0,5) = -Y(u,v = 0,r = 0,-5) (50) 

Hence, in analogy with the previous section, the even derivatives of Y with respect to v (and 
then r, and 5) must be zero, or: 



Y 


= 0; 


Y 

1 vvu 


= 


Y 


= 0; 


Y 

rru 


= 


Y ss 


= 0; 


Y&6u 


= 



(51) 



Finally, due to port/starboard symmetry, the force Y should not be affected by u when 
v = 0, r = 0, 5 = 0, except for propeller effects, which break the symmetry. For this reason, 
we will include terms such as Y u to allow for the propeller asymmetry (twin propeller ships 
with counter-rotating propellers have zero such terms). 
Finally, we derive the following expansion for Y: 



Y = Y e + Y u u + Y uu u 2 + Y v v + Y r r + Y v v + Y r r + Y s 5 + Y Su 5u + Y vu vu + 

Y rr5 r 2 5 + Y rrv r 2 v + Y vvr v 2 r + Y vvS v 2 5 + Y vr $vr5 + Y S s r 5 2 r + Y$ Sv 5 2 v (52) 
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3.2.3 Fluid Moment N 

The derivation for N follows the same exact steps as for the side force Y, i.e. the same 
symmetries apply. As a result, we first find that the following coefficients are zero: 



N r 



N ss = 



N, 



= 

, . a 



N SS u = 



(53) 



Hence, we derive an analogous expansion for the moment: 



N = N e + N u u + N uu u 2 + N^v + Nf.f + N v v + N r r + N S S + N Su Su + 

N vu vu + N ru ru + N vuu vu 2 + N ruu ru 2 + N 5uu 5u 2 + N vvv v 3 + A^ rrr r 3 + 

NsssS 3 + N rr5 r 2 S + N rrv r 2 v + N vvr v 2 r + N vv5 v 2 S + N vrS vrS + 

Ns Sr S 2 r + NssvPv. (54) 



4 VESSEL DYNAMICS: LINEAR CASE 

4.1 Surface Vessel Linear Model 

We first discuss some of the hydrodynamic parameters which govern a ship maneuvering in 
the horizontal plane. The body x-axis is forward and the y-axis is to port, so positive r has 
the vessel turning left. We will consider motions only in the horizontal plane, which means 
9 = ip = p = q = w = 0. Since the vessel is symmetric about the x — z plane, yc = 0; zq is 
inconsequential. We then have at the outset 



(du 2 \ 
— — rv — xqt I (55) 

/ dv dr\ 
Y = m^ + ru + xa-j 

dr f dv \ 

= I -m +mXG {m +ru )- 

Letting u — U + u, where U >> u, and eliminating higher-order terms, this set is 



X = m | (56) 
( 9v TT dr\ 

dr f dv \ 



18 



4 VESSEL DYNAMICS: LINEAR CASE 



A number of coefficients can be discounted, as noted in the last chapter. First, in a homoge- 
neous sea, with no current, wave, or wind effects, {X x , X y , X^, Y x , Y y , Y^, N x , N y , N^} are all 
zero. We assume that no hydrodynamic forces depend on the position of the vessel. 1 Second, 
consider X v : since this longitudinal force would have the same sign regardless of the sign of 
v (because of side-to-side hull symmetry), it must have zero slope with v at the origin. Thus 
X v = 0. The same argument shows that {X r , Xy, X r , Y u , N u , N^} = 0. Finally, since 
fluid particle acceleration relates linearly with pressure or force, we do not consider nonlin- 
ear acceleration terms, or higher time derivatives. It should be noted that some nonlinear 
terms related to those we have eliminated above are not zero. For instance, Y uu = because 
of hull symmetry, but in general X vv = only if the vessel is bow-stern symmetric. 
We have so far, considering only the linear hydrodynamic terms, 



(m-Xu)u = X u u + X' (57) 
(m - Y v )v + (mx G - Yr)r = Y v v + (Y r — mU)r + Y' (58) 
(mx G - N v )i) + (I zz - N r )r = N v v - (N r - mx G U)r + N'. (59) 

The right side here carries also the imposed forces from a thruster(s) and rudder(s) {X', Y', N'} 
Note that the surge equation is decoupled from the sway and yaw, but that sway and yaw 
themselves are coupled, and therefore are of immediate interest. With the state vector 
s — {v,r} and external force/moment vector F = {Y',N'}, a state-space representation of 
the sway/yaw system is 



m — Yy 
mx G — Ni 



mx G — Y^ 



ds Y v Y r - mU 

dt N v N r - mx G U 

M's = Ps + F 

's = M~ 1 Ps + M~ l F 

's = As + BF. 



s + F, or 



(60) 



(61) 



The matrix M is a mass or inertia matrix, which is always invertible. The last form of the 
equation is a standard one wherein A represents the internal dynamics of the system, and 
B is a gain matrix for the control and disturbance inputs. 



4.2 Stability of the Sway /Yaw System 

Consider the homogeneous system s — As: 



si = A ll s l + A l2 s 2 

S 2 — A21S1 + ^22^2- 

1 Note that the linearized heave/pitch dynamics of a submarine do depend on the pitch angle; this topic 
will be discussed later. 



4.2 Stability of the Sway/Yaw System 
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We can rewrite the second equation as 



'do 

dt 

and substitute into the first equation to give 



-l 



S2= [dt~ A22 ) MlSl (62) 



si + (-An - ^ 22 )si + (AhAm - A 12 A 21 ) Sl = 0. (63) 

Note that these operations are allowed because the derivative operator is linear; in the lan- 
guage of the Laplace transform, we would simply use s. A necessary and sufficient condition 
for stability of this ODE system is that each coefficient must be greater than zero: 



-A u -A 22 > (64) 
A U A 22 - A l2 A 2l > 



The components of A for the sway/yaw problem are 



(I zz - N,)Y V + (Yr - mx G )N v 
{m-Yi){I zz - Nr)-{mx G -Yr){mx G - Ni) 

-{I zz - N r )(mU - Y r ) - (Yr - mx G )(mx G U - N r ) 
(m - Yy)(I zz - N r ) - (mx G - Y r )(mx G - Ny) 

(N v - mx G )Y v + (m - Y V )N V 

(m - Y V )(I ZZ - N r ) - (mx G - Y r )(mx G - Ny) 

-(Nj, - mx G )(mU - Y r ) - (m - Yt)(mx G U - N r) 
(m - Yv)(I zz - N r ) - (mx G - Y r )(mx G - N*) 

The denominators are identical, and can be simplified. First, let x G ~ 0; valid for many 
vessels with the origin is at the geometric center. If the origin is at the center of mass, 
x G = 0. Next, if the vessel is reasonably balanced with regard to forward and aft areas with 
respect to the origin, the terms {Ny,Y r ,N v ,Y r } take very small values in comparison with 
the others. To wit, the added mass term — Y^ is of the order of the vessel's material mass 
m, and similarly N r ~ —I zz - Both Y% and N r take large negative values. Linear drag and 
rotational drag are significant also; these are the terms Y v and N r , both large and negative. 
The denominator for A's components reduces to (m — Yu)(I zz — N r ), and 



A n = 

A 12 = 

A 21 = 

A 22 = 



^ = £yr < 

A 22 = ^ < o. 

Hence the first condition for stability is met: —An — A 22 > 0. For the second condition, 
since the denominators of the are identical, we have only to look at the numerators. For 
stability, we require 
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(I^-N+Wm-YjNr (66) 
- [NyY v + {m- Y V )N V ] [-{I zz - N,){mU - Y r ) + Y,N r ] > 0. 

The first term is the product of two large negative and two large positive numbers. The 
second part of the second term contains mil, which has a large positive value, generally 
making stability critical on the (usually negative) N v . When only the largest terms are 
considered for a vessel, a simpler form is common: 

C = Y v N r + N v (mU-Y r )>0. (67) 

C is called the vessels stability parameter. The terms of C compete, and yaw/sway stability 
depends closely on the magnitude and sign of N v . Adding more surface area aft drives N v 
more positive, increasing stability as expected. Stability can also be improved by moving 
the center of gravity forward. Nonzero x G shows up as follows: 

C = Y v (N r - mx G U) + N v (mU - Y r ) > 0. (68) 

Since iV r and Y v are both negative, positive x G increases the (positive) influence of C"s first 
term. 




Figure 3: Convention for positive rudder angle in the vessel reference system. 



4.3 Basic Rudder Action in the Sway /Yaw Model 

Rudders are devices which develop large lift forces due to an angle of attack with respect 
to the oncoming fluid. The basic form is as follows: L = ^pAU 2 Ci(a), where a is the angle 
of attack. The lift coefficient Ci is normally linear with a near a = 0, but the rudder stalls 
when the angle of attack reaches a critical value, and thereafter develops much less lift. We 
will assume that a is small enough that the linear relationship applies: 



Ci(a) 



da 
mo: 

distance I aft, the moment equation is quite simple. We have 



a. (69) 



=o 

Since the rudder develops force (and a small moment) far away from the body origin, say a 



4.3 Basic Rudder Action in the Sway/Yaw Model 
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IU 2 . (71) 

Note the difference between the rudder angle expressed in the body frame, 5, and the total 
angle of attack a. Angle of attack is influenced by <5, as well as v/U and Ir. Thus, in tank 
testing with v = 0, S = a and N$ = N a , etc., but in real conditions, other hydrodynamic 
derivatives are augmented to capture the necessary effects, for example N v and N r . Generally 
speaking, the hydrodynamic characteristics of the vessel depend strongly on the rudder, even 
when 5 = 0. In this case the rudder still opposes yaw and sway perturbations and acts to 
stabilize the vessel. 

A positive rudder deflection (defined to have the same sense as the yaw angle) causes a 
negative yaw perturbation, and a very small positive sway perturbation. 

4.3.1 Adding Yaw Damping through Feedback 

The stability coefficient C resulting from the addition of a control law 5 = k r r, where k r > 
is a feedback gain, is 

C = Y v (N r - mx G U + k r N s ) + N v (mU -Y r - k r Y 5 ). (72) 

Ys is small positive, but N$ is large and negative. Hence C becomes more positive, since Y v 
is negative. 

Control system limitations and the stalling of rudders make obvious the fact that even a very 
large control gain k r cannot completely solve stability problems of a poorly-designed vessel 
with an inadequate rudder. On the other hand, a vessel which is overly stable (C >> with 
no rudder action) is unmaneuverable. A properly-balanced vessel just achieves stability with 
zero rudder action, so that a reasonable amount of control will provide good maneuvering 
capabilities. 

4.3.2 Heading Control in the Sway/Yaw Model 

Considering just the yaw equation of motion, i.e., v — 0, with a rudder, we have 

(I zz - Nf)(f> + (mx G U - N r )j) = N S S. (73) 

Employing the control law 5 = k^cfi, the system equation becomes a homogeneous, second- 
order ODE: 

(I zz - Nr)4> + (mx G U - N r )<j) - Ngk^ = 0. (74) 

Since all coefficients are positive (recall N$ < 0), the equation gives a stable 9 response, 
settling under second-order dynamics to 9{oo) = 0. The control law 5 = k^{(j) — (fidesired) + k r r 
is the basis for heading autopilots, which are used to track (fidesired- This use of an error signal 



U 2 



(70) 
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to drive an actuator is in fact the essence of feedback control. In this case, we require sensors 
to obtain r and 0, a controller to calculate 5, and an actuator to implement the corrective 
action. 

4.4 Response of the Vessel to Step Rudder Input 
4.4.1 Phase 1: Accelerations Dominate 

When the rudder first moves, acceleration terms dominate, since the velocities are zero. The 
equation looks like this: 



m — Yi, mxc — Yf 
mx G - Ni I zz - iV, 



r 



M = \ x 5 

r 1 N s 



S. (75) 



Since Y+ and mxc are comparatively small in the first row, we have 



*< > = A < 76 > 



and the vessel moves to the left, the positive ^-direction. The initial yaw is in the negative 
r-direction, since Ng < 0: 

± zz ly r 

The first phase is followed by a period (Phase 2), in which many terms are competing and 
contributing to the transient response. 

4.4.2 Phase 3: Steady State 

When the transients have decayed, the vessel is in a steady turning condition, and the 
accelerations are zero. The system equations reduce to 



(78) 



J v \ - - I ( mx cU - N r )Y s + (Y r - mU)N s \ 
\r )- C { N V Y S - Y V N S ) ' 

Note that the denominator is the stability parameter. The steady turning rate is thus 
approximated by 



YvNs -5. (79) 



C 

With C > 0, the steady-state yaw rate is negative. If the vessel is unstable (C < 0), it 
turns in the opposite direction than expected. This turning rate equation can also be used 
to estimate turning radius R: 

r, U UC 

The radius goes up directly with C, indicating that too stable a ship has poor turning 
performance. We see also that increasing the rudder area increases N$, decreasing R as 
desired. Increasing the deflection 5 to reduce R works only to the point of stalling. 



4.5 Summary of the Linear Maneuvering Model 
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4.5 Summary of the Linear Maneuvering Model 

We conclude our discussion of the yaw/sway model by noting that 

1. The linearized sway/yaw dynamics of a surface vessel are strongly coupled, and they 
are independent of the longitudinal dynamics. 

2. The design parameter C should be slightly greater than zero for easy turning, and " 
hands-off" stability. The case C < should only be considered under active feedback 
control. 

3. The analysis is valid only up to small angles of attack and turning rates. Very tight ma- 
neuvering requires the nonlinear inertial components and hydrodynamic terms. Among 
other effects, the nonlinear equations couple surge to the other motions, and the actual 
vessel loses forward speed during maneuvering. 

4.6 Stability in the Vertical Plane 

Stability in the horizontal plane changes very little as a function of speed, because drag and 
lift effects generally scale with U 2 . This fact is not true in the vertical plane, for which the 
dimensional weight/buoyancy forces and moments are invariant with speed. For example, 
consider the case of heave and pitch, with xq = and no actuation: 



Z^w + Z w w + Z q q + Z q q + (B - W) (81) 
M^w + M w w + M q q + M q q - Bl b sin 9. (82) 

The last term in each equation is a hydrostatic effect induced by opposing net buoyancy 
B and weight W . l b denotes the vertical separation of the center of gravity and the center 
of buoyancy, creating the so-called righting moment which nearly all underwater vehicles 
possess. Because buoyancy effects do not change with speed, the dynamic properties and 
hence stability of the vehicle may change with speed. 

5 SIMILITUDE 

5.1 Use of Nondimensional Groups 

For a consistent description of physical processes, we require that all terms in an equation 
must have the same units. On the basis of physical laws, some quantities are dependent on 
other, independent quantities. We form nondimensional groups out of the dimensional ones 
in this section, and apply the technique to maneuvering. 

The Buckingham 7r-theorem provides a basis for all nondimensionalization. Let a quantity 
Q n be given as a function of a set of n — 1 other quantities: 



m 



Uq 



dq 
yy dt 
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Qn — /q(<5i, <?2, * • * , Qn-l)- (83) 

There are n variables here, but suppose also that there are only k independent ones; k is 
equivalent to the number of physical unit types encountered. The theorem asserts that there 
are n — k dimensionless groups 7Tj that can be formed, and the functional equivalence is 
reduced to 

^n-k — U(^l,^2,'",^n-k-l)- (84) 

Example. Suppose we have a block of mass m resting on a frictionless horizontal surface. 
At time zero, a steady force of magnitude F is applied. We want to know X(T), the distance 
that the block has moved as of time T. The dimensional function is X(T) = /q(jto, F, T), so 
n — 4. The (MKS) units are 

[*(•)] = rn 
[m] = kg 

[F] = kgm/s 2 
[T] = s, 

and therefore k — 3. There is just one nondimensional group in this relationship; m assumes 
only a constant (but unknown) value. Simple term-cancellation gives 7Ti = X(T)m/FT 2 , 
not far at all from the known result that X(T) = FT 2 /2m\ 

Example. Consider the flow rate Q of water from an open bucket of height h, through a 
drain nozzle of diameter d. We have 

Q = fQ(h,d,p,fi,g), 

where the water density is p, and its absolute viscosity p; g is the acceleration due to gravity. 
No other parameters affect the flow rate. We have n — 6, and the (MKS) units of these 
quantities are: 

[Q] = m 3 /s 
[h] = m 
[d] = m 
[p] = kg/m 3 
[p] = kg/ms 
[g] = m/s 2 

There are only three units that appear: [length, time, mass], and thus k — 3. Hence, 
only three non-dimensional groups exist, and one is a unique function of the other two. To 
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arrive at a set of valid groups, we must create three nondimensional quantities, making sure 
that each of the original (dimensional) quantities is represented. Intuition and additional 
manipulations come in handy, as we now show. 

Three plausible first groups are: 7Ti = pQ/dfj,, 7r 2 = dp^/gh/p, and 7r 3 = h/d. Note that all 
six quantities appear at least once. Since h and d have the same units, they could easily 
change places in the first two groups. However, ix\ is recognized as a Reynolds number 
pertaining to the orifice flow. 7T2 is more awkward, but products and fractions of groups are 
themselves valid groups, and we may construct — t^\/t^2 — Q/d 2 ^/gh to nondimensionalize 
Q with a pressure velocity, and then 7T 5 = tti/tt± — pd^fghj '[/, to establish an orifice Reynolds 
number independent of Q. We finally have the useful result 



7T4 = fn(n 5 , 7T 2 ) > 

Q = f ( P d ^9h h\ 
d 2 y/gh n \ /i ' d J 

The uncertainty about where to use h and d, and the questionable importance of h/d as a 
group are remnants of the theorem. Intuition is that h/d is immaterial, and the other two 
terms have a nice physical meaning, e.g., 7r 5 is a Reynolds number. 

The power of the 7r-theorem is primarily in reducing the number of parameters which must 
be considered independently to characterize a process. In the flow example, the theorem 
reduced the number of independent parameters from five to two, with no constraints about 
the actual physics taking place. 

5.2 Common Groups in Marine Engineering 

One frequently encounters the following groups in fluid mechanics and marine engineering: 
1. Froude number: 



where U is the speed of the vessel, g is the acceleration due to gravity, and L is 
the waterline length of the vessel. The Froude number appears in problems involving 
pressure boundary conditions, such as in waves on the ocean surface. Roughly speaking, 
it relates the vessel speed U to water wave speeds of wavelength L; the phase speed of 
a surface wave is 7r, where A is the wavelength. 

2. Cavitation number: 

where represents the ambient total pressure, P v the vapor pressure of the fluid, 
and U the propeller inlet velocity. A low cavitation number means that the Bernoulli 
pressure loss across the lifting surface will cause the fluid to vaporize, causing bubbles, 
degradation of performance, and possible deterioration of the material. 
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3. Reynolds number: 



Re 



Ul 
pjp 



(87) 



where U is velocity, p is absolute viscosity, and p is density. Since Re appears in many 
applications, I represents one of many length scales. Reynolds number is a ratio of fluid 
inertial pressures to viscous pressures: When Re is high, viscous effects are negligible. 
Re can be used to characterize pipe flow, bluff body wakes, and flow across a plate, 
among others. 

4. Weber number: 



W 



pU 2 l 



where a is the surface tension of a fluid. Given that [a] = N/m (MKS), pU 2 normalizes 
pressure, and I normalizes length. The Weber number indicates the importance of 
surface tension. 

To appreciate the origins of these terms from a fluid particle's point of view, consider a box 
having side lengths [dx, dy, dz\. Various forces on the box scale as 



pgdxdydz ~ pgl 3 
Pdxdy ~ PI 2 



(inertia) F; 

(gravity) F g 
(pressure) F p 

(shear) F s 
(surface tension) F a = adx ~ al. 

Thus the groups listed above can be written as 



dv 7 7 7 9v . . . tt2 .o 
p—dxdyaz + pv—dxdydz ~ pu I 
ot ox 



p^-dxdy 

oz 



pXJl 



— ~ 



Re = 



W 
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u 2 




gi 
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c v 
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pU 2 


F 


pUl 


~F 

1 s 
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F 


P U 2 l 


± a 


a 



When testing models, it is imperative to maintain as many of the nondimensional groups as 
possible of the full-scale system. This holds for the geometry of the body and the kinematics 
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of the flow, the surface roughness, and all of the relevant groups governing fluid dynamics. 
Consider the example of nozzle flow from a bucket. Suppose that we conduct a model 
test in which Re is abnormally large, i.e., the viscous effects are negligible. Under inviscid 
conditions, the flow rate is Q = ird 2 y/2gh/ '4. This rate cannot be achieved for lower- Re 
conditions because of fluid drag in the orifice, however. 

In a vessel, we write the functional relationship for drag as a starting point: 

C r = = fQ(P>P>9,<r,U,l) 

= U(Re,Fr,W). 

First, since I is large, W is very large, and hence surface tension plays no role. Next, we look 
at Re — Ul/v and Fr = U/y/gl, both of which are important for surface vessels. Suppose 
that l s hip = XI model, so that usually A >> 1; additionally, we set g mo dei = 9 ship, i-e., the model 
and the true vessel operate in the same gravity field. 

Froude number similitude requires U mo d e i = U s hi P /VX. Then Reynolds number scaling im- 
plies directly v mode i = v shiv IX z l 2 . Unfortunately, few fluids with this property are workable 
in a large testing tank. As a result, accurate scaling of Re for large vessels to model scale is 
quite difficult. 

For surface vessels, and submarines near the surface, it is a routine procedure to employ 
turbulence stimulators to achieve flow that would normally occur with ship-scale Re. Above 
a critical value Re ~ 500, 000, C/ is not sensitive to Re. With this achieved, one then tries 
to match Fr closely. 

5.3 Similitude in Maneuvering 

The linear equations of motion for the horizontal yaw/sway problem are: 

(m - Yi)v - Y v v + (mil - Y r )r + (mx G - Y f )r = Y 
(I zz - Nf.)r + (mx G U - N r )r + (mx G - N*)v - N v v = N. 



These equations can be nondimensionalized in a standard way, by using the quantities 
[U,L,p]: these three values provide the necessary units of length, time, and mass, and fur- 
thermore are readily accessible to the user. First, we create nondimensional states, denoted 
with a prime symbol: 



v' = 
v> = 

¥ = 



(89) 



28 



5 SIMILITUDE 



L 

r = _ r . 

We follow a similar procedure for the constant terms as follows, including a factor of 1/2 
with p, for consistency with our previous expressions: 

m! = y~T^ (9°) 



I' 

zz 



Xq 



U' 
Y- 

Y' 

Y! 

± r 

Y' 
Y' 
N' 
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N' 
N' 



\ P U> 

Izz 

w 

L 

U 

Y v 
\pUL 2 

Y r 



\pUU> 
Y 



\pU 2 L 2 



N v 



\pUL* 

Nr 



N r 

jpur* 

N 



Note that every force has been normalized with ^pU 2 L 2 , and every moment with ^pU 2 L 3 ; 
time has been also nondimensionalized with L/U. Thus we arrive at a completely equivalent 
set of nondimensional system equations, 

(rri - Y!)v' - y>' + (m'U' - Y' r )r' + (m'x' G - Y!)r' = Y' (91) 
(I' zz - N^r' + (m'x' G U' - N' r )r' + (m'x' G - N<)v' - N' v v' = N'. 
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Since fluid forces and moments generally scale with U 2 , the nondimensionalized description 
holds for a range of velocities. 

5.4 Roll Equation Similitude 

Certain nondimensional coefficients may arise which depend explicitly on U, and therefore 
require special attention. Let us carry out a similar normalization of the simplified roll 
equation 

(I xx - K p )p - K pP - K^> = K. (92) 

For a surface vessel, the roll moment is based on metacentric stability, and has the form 
K^j = —pgV(GM), where V is the displaced fluid volume of the vessel, and GM is the 
metacentric height. The nondimensional terms are 

r - = w (93) 



K 



p' = TTP 



K p 
\pU> 
K p 
\pUL* 

±pU 2 L 3 

u 2 

L 



leading to the equivalent system 

(I'xx ~ K)P' ~ KP' ~ (^) V'(GM')^ = K>. (94) 

Note that the roll angle (p was not nondimensionalized. The Froude number has a very strong 
influence on roll stability, since it appears explicitly in the nondimensional righting moment 
term, and also has a strong influence on K' p . 

In the case of a submarine, the righting moment has the form = —Bh, where B is the 
buoyant force, and h is the righting arm. The nondimensional coefficient becomes 



K> = — 
^ \p\J 2 I?' 

K[p again depends strongly on U, since B and h are fixed; this K'^ needs to be maintained 
in model tests. 
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6 CAPTIVE MEASUREMENTS 

Before making the decision to measure hydrodynamic derivatives, a preliminary search of 
the literature may turn up useful estimates. For example, test results for many hull-forms 
have already been published, and the basic lifting surface models are not difficult. The 
available computational approaches should be considered as well; these are very good for 
predicting added mass in particular. Finally, modern sensors and computer control systems 
make possible the estimation of certain coefficients based on open-water tests of a model or 
a full-scale design. 

In model tests, the Froude number Fr = which scales the influence of surface waves, 
must be maintained between model and full-scale surface vessels. Reynolds number Re = 
which scales the effect of viscosity, need not be matched as long as the scale model 
attains turbulent flow (supercritical Re). One can use turbulence stimulators near the bow 
if necessary. Since the control surface(s) and propeller(s) affect the coefficients, they should 
both be implemented in model testing. 

6.1 Towtank 

In a towtank, tow the vehicle at different angles of attack, measuring sway force and yaw 
moment. The slope of the curve at zero angle determines Y v and N v respectively; higher- 
order terms can be generated if the points deviate from a straight line. Rudder derivatives 
can be computed also by towing with various rudder angles. 

6.2 Rotating Arm Device 

On a rotating arm device, the vessel is fixed on an arm of length R, rotating at constant 
rate r: the vessel forward speed is U = rR. The idea is to measure the crossbody force and 
yaw moment as a function of r, giving the coefficients Y r and N r . Note that the lateral force 
also contains the component (m — Yy)r 2 R. The coefficients Y v and N v can also be obtained 
by running with a fixed angle of attack. Finally, the measurement is made over one rotation 
only, so that the vessel does not re-enter its own wake. 

6.3 Planar-Motion Mechanism 

With a planar motion mechanism, the vessel is towed at constant forward speed U, but 
is held by two posts, one forward and one aft, which can each impose independent sway 
motions, therefore producing variable yaw. The model moves in pure sway if y a {t) = Vb(t), 
in a pure yaw motion about the mid- length point if y a (t) = —yb(t), or in a combination sway 
and yaw motion. The connection points are a distance I forward and aft from the vessel 
origin. 

Usually a sinusoidal motion is imposed: 



y a (t) = a cos cut 
y b (t) = bcos(cot + i/j), 



(95) 
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and the transverse forces on the posts are measured and approximated as 



Y a (t) = F a cos(ut + 6 a ) (96) 
Y b (t) = F b cos(ut + 6 b ). 



If linearity holds, then 



(m - Yy)v - Y v v + (mU - Y r )r + (mx G - Y+)r = Y a + Y b (97) 
(I zz - Nr)r + (mx G U - N r )r + {mx G - N^v - N v v = (Y b -Y a )l. 

We have v — (y a + y a )/2 and r = (y b — y a )/^l- When a = b, these become 

v = — ^ (sino;t(l + cos-0) + coswt sin-0) (98) 

au 2 . . 

v = {cosoot{l + cosip) — smut smip) 

2 

CLUJ . . . . ■ i \ 

r = (smut^osip — 1) + cos out smip) 

Zii 

au 2 . . ... 

r = — [cos uot{cos ip — 1) — sinwt smip) . 

Equating the sine terms and then the cosine terms, we obtain four independent equations: 

(m-Y,)(-^(l + cost/0- (99) 
Y v i — —J smip + 



(mU-Y r ) (-^) sin^ + 



(mx G - Yf.) I — — I (cosV> - 1) = F a cos 9 a + F b cos 9 b 



(m-Y v )[-^Pj(-smip)- (100) 

Y v 2"J (1 + cost/0 + 
(mU-Y r ) (-^) (oos^-l) + 

(mx G - Yr) (y- (- sin ip) = -F a sin 8 a - F b sin 9 b 
(Izz-N*) (~^P) (cosV-l)+ (101) 



32 



6 CAPTIVE MEASUREMENTS 



(mx G U - N r ) (-7^-) sinV' + 
(mx G - Ny) ( - C ^— ) (1 + cos if)) - 



N v (--y) sin V' = K-Pft cos - COS 6 a ) 



{I zz -N,)\-—U-s^)+ (102) 



(mx G t7 - N r ) (-7^) (cos^ - 1) + 



aw 2 ' 



(m^G - I -— 1 (-sin-0) - 



N v (- y) t 1 + cos ^) = l (~ Fh sin ^ + Fa sin 9 *) 



In this set of four equations, we know from the imposed motion the values [U, if), a, uj\. From 
the experiment, we obtain [F a , Fb,0 a , 8b], and from the rigid-body model we have [m, I zz , xq\. 
It turns out that the two cases of ip = (pure sway motion) and ip = 180° (pure yaw 
motion) yield a total of eight independent equations, exactly what is required to find the 
eight coefficients [Y v , Y v , Y r , Y r , N v , N v , N+, N r ]. Remarkably, we can write the eight solutions 
directly: For ip = 0, 



au) 2 ^ 



(m - Y-o) I —I (2) = F a cos 9 a + F b cos 6 b (103) 



-Y v (-™^(2) = -F a sm9 a -F b sm9 b 
(mx G - Nt) (-^r) (2) = l(F b cos 9 b - F a cos6 a ) 



-N v (-™^(2) = l(-F h smd b + F a smd a ), 
to be solved respectively for [Y$, Y v ,Ni, N v ]. For -0 = 180°, we have 



(mx G — Yj.) I — — j (—2) = F a cos 9 a + F b cos 9 b (104) 



(mC/-y r )(-^)(-2) = -F a sm9 a ~F b sin9 b 



au) 2 ^ 



(I zz -Nr)^-—j(-2) = l(F b cos 9 b -F a cos 9 a ) 
(mx G U-N r ) (-2) = J(-F b sin0 6 + F a sin0 a ), 
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to be solved for [Y+, Y r ,Nr, N r }. Thus, the eight linear coefficients for a surface vessel maneu- 
vering, for a given speed, can be deduced from two tests with a planar motion mechanism. 
We note that the nonlinear terms will play a significant role if the motions are too large, 
and that some curve fitting will be needed in any event. The PMM can be driven with more 
complex trajectories which will target specific nonlinear terms. 

7 STANDARD MANEUVERING TESTS 

This section describes some of the typical maneuvering tests which are performed on full-scale 
vessels, to assess stability and performance. 

7.1 Dieudonne Spiral 

1. Achieve steady speed and direction for one minute. No changes in speed setting are 
made after this point. 

2. Turn rudder quickly by 15°, and keep it there until steady yaw rate is maintained for 
one minute. 

3. Reduce rudder angle by 5°, and keep it there until steady yaw rate is maintained for 
one minute. 

4. Repeat in decrements of -5°, to -15°. 

5. Proceed back up to 15°. 

The Dieudonne maneuver has the vessel path following a growing spiral, and then a con- 
tracting spiral in the opposite direction. The test reveals if the vessel has a memory effect, 
manifested as a hysteresis in yaw rate r. For example, suppose that the first 15° rudder 
deflection causes the vessel to turn right, but that the yaw rate at zero rudder, on the way 
negative, is still to the right. The vessel has gotten "stuck" here, and will require a negative 
rudder action to pull out of the turn. But if the corrective action causes the vessel to turn left 
at all, the same memory effect may occur. It is easy to see that the rudder in this case has 
to be used excessively driving the vessel back and forth. We say that the vessel is unstable, 
and clearly a poor design. 

7.2 Zig-Zag Maneuver 

1. With zero rudder, achieve steady speed for one minute. 

2. Deflect the rudder to 20°, and hold until the vessel turns 20°. 

3. Deflect the rudder to -20°, and hold until the vessel turns to -20o with respect to the 
starting heading. 

4. Repeat. 
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This maneuver establishes several important characteristics of the yaw response. These are: 
the response time (time to reach a given heading), the yaw overshoot (amount the vessel 
exceeds ±20° when the rudder has turned the other way), and the total period for the 20° 
oscillations. Of course, similar tests can be made with different rudder angles and different 
threshold vessel headings. 

7.3 Circle Maneuver 

From a steady speed, zero yaw rate condition, the rudder is moved to a new setting. The 
vessel responds by turning in a circle. After steady state is reached again, parameters of 
interest are the turning diameter, the drift angle (3, the speed loss, and the angle of heel ip. 

7.3.1 Drift Angle 

The drift angle is the equivalent to angle of attack for lifting surfaces, and describes how the 
vessel "skids" during a turn. If the turning circle has radius R (measured from the vessel 
origin), then the speed tangential to the circle is U — rR. The vessel-reference velocity 
components are thus u = U cos/3 and v — — Usia/3. A line along the vessel centerline 
reaches closest to the true center of the turning circle at a point termed the turning center. 
At this location, which may or may not exist on the physical vessel, there is no apparent 
lateral velocity, and it appears to an observer there that the vessel turns about this point. 

7.3.2 Speed Loss 

Speed loss occurs primarily because of drag induced by the drift angle. A vessel which drifts 
very little may have very little speed loss. 

7.3.3 Heel Angle 

Heel during turning occurs as a result of the intrinsic coupling of sway, yaw, and roll caused 
by the center of gravity. In a surface vessel, the fluid forces act below the waterline, but the 
center of gravity is near the waterline or above. When the rudder is first deflected, inertial 
terms dominate (Phase 1) and the sway equation is 

(m - Y^v - (Yr - mx G )r = Y s 5. (105) 

The coefficients for r are quite small, and thus the vessel first rolls to starboard (positive) 
for a positive rudder action. 

When steady turning conditions are reached (Phase 3), hydrodynamic forces equalize the 
centrifugal force mUr and the rudder force Yg8. The sway equation is 

-Y v v + (mU - Y r )r = Y S S, (106) 

with Y r small, v < when r > for most vessels, and \Y v v\ > \Ys6\. Because the centrifugal 
force acts above the waterline, the vessel ultimately rolls to port (negative) for a positive 
rudder action. 
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The transition between the inertially-dominated and steady-turning regimes includes an 
overshoot: in the above formulas, the vessel overshoots the final port roll angle as the vessel 
slows. From the sway equation, we see that if the rudder is straightened out at this point, 
the roll will momentarily be even worse! 

In summary, the vessel rolls into the turn initially, but then out of the turn in the steady 
state. 

7.3.4 Heeling in Submarines with Sails 

Submarines typically roll into a turn during all phases. Unlike surface vessels, which have 
the rigid mass center above the center of fluid forcing, submarines have the mass center 
below the rudder action point, and additionally feel the effects of a large sail above both. 
The inert ial equation 

(m - Y v )v - (Yf - mx G )r = Y s 5 (107) 

is dominated by mv (acting low), —YyV (acting high), and Y$5 (intermediate). Because 
\Y V \ > m, the vessel rolls under the sail, the keel out of the turn. In the steady state, 

-Y v v + [mU - Y r )r = Y s 5. (108) 

The drift angle (3 keeps the Yy-force, acting high, toward the center of the turn, and again 
centrifugal force mUr causes the bottom of the submarine to move out of the turn. Hence, 
the roll angle of a submarine with a sail is always into the turn, both initially and in the 
steady state. The heel angle declines as the speed of the submarine drops. 

8 STREAMLINED BODIES 

8.1 Nominal Drag Force 

A symmetric streamlined body at zero angle of attack experiences only a drag force, which 
has the form 

F A = - 1 - P C A A U 2 . (109) 

The drag coefficient C A has both pressure and skin friction components, and hence area A Q 
is usually that of the wetted surface. Note that the A-subscript will be used to denote zero 
angle of attack conditions; also, the sign of F A is negative, because it opposes the vehicle's 
x-axis. 

8.2 Munk Moment 

Any shape other than a sphere generates a moment when inclined in an inviscid flow. 
d'Alembert's paradox predicts zero net force, but not necessarily a zero moment. This 
Munk moment arises for a simple reason, the asymmetric location of the stagnation points, 
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where pressure is highest on the front of the body (decelerating flow) and lowest on the back 
(accelerating flow). The Munk moment is always destabilizing, in the sense that it acts to 
turn the vehicle perpendicular to the flow. 

Consider a symmetric body with added mass components A xx along the vehicle (slender) x- 
axis (forward), and A zz along the vehicle's z-axis z (up). We will limit the present discussion 
to the vertical plane, but similar arguments can be used to describe the horizontal plane. 
Let a represent the angle of attack, taken to be positive with the nose up - this equates 
to a negative pitch angle (p in vehicle coordinates, if it is moving horizontally. The Munk 
moment is: 



M m = - - (A zz - A XX )U 2 sin 2a (110) 
— ~(A ZZ — A xx )U 2 a. 

A zz > A xx for a slender body, and the negative sign indicates a negative pitch with respect to 
the vehicle's pitch axis. The added mass terms A zz and A xx can be estimated from analytical 
expressions (available only for regular shapes such as ellipsoids), from numerical calculation, 
or from slender body approximation (to follow). 



8.3 Separation Moment 

In a viscous fluid, flow over a streamlined body is similar to that of potential flow, with the 
exceptions of the boundary layer, and a small region near the trailing end. In this latter 
area, a helical vortex may form and convect downstream. Since vortices correlate with low 
pressure, the effect of such a vortex is stabilizing, but it also induces drag. The formation 
of the vortex depends on the angle of attack, and it may cover a larger area (increasing the 
stabilizing moment and drag) for a larger angle of attack. For a small angle of attack, the 
transverse force F v can be written in the same form as for control surfaces: 



F n = \ P C n A U 2 (111) 
K 2^ aA ' U - 

With a positive angle of attack, this force is in the positive ^-direction. The zero-a drag 
force Fa is modified by the vortex shedding: 



~ P C a A U 2 , where (112) 



C a = C A cos 2 cf). 



The last relation is based on writing Ca(U cos^) 2 as (Ca cos 2 (f))U 2 , i.e., a decomposition 
using apparent velocity. 
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8.4 Net Effects: Aerodynamic Center 

The Munk moment and the moment induced by separation are competing, and their mag- 
nitudes determine the stability of a hullform. First we simplify: 

Fa 7a 
F n = 7n« 

M m = -7 m a. 

Each constant 7 is taken as positive, and the signs reflect orientation in the vehicle reference 
frame, with a nose-up angle of attack. The Munk moment is a pure couple which does not 
depend on a reference point. We pick a temporary origin O for F n however, and write the 
total pitch moment about O as: 

M = M m + FJ n (113) 

= (~7m + JnQa. 

where l n denotes the (positive) distance between O and the application point of F n . The net 
moment about O is zero if we select 

L = ^, (H4) 

In 

and the location of O is then called the aerodynamic center or AC. 

The point AC has an intuitive explanation: it is the location on the hull where F n would act 
to create the total moment. Hence, if the vehicle's origin lies in front of AC, the net moment 
is stabilizing. If the origin lies behind AC, the moment is destabilizing. For self-propelled 
vehicles, the mass center must be forward of AC for stability. Similarly, for towed vehicles, 
the towpoint must be located forward of AC. In many cases with very streamlined bodies, 
the aerodynamic center is significantly ahead of the nose, and in this case, a rigid sting would 
have to extend at least to AC in order for stable towing. As a final note, since the Munk 
moment persists even in inviscid flow, AC moves infinitely far forward as viscosity effects 
diminish. 

8.5 Role of Fins in Moving the Aerodynamic Center 

Control surfaces or fixed fins are often attached to the stern of a slender vehicle to enhance 
directional stability. Fixed surfaces induce lift and drag on the body: 

L=^pA f U 2 d(a) ~ lL a (115) 
D = —^-pAfU 2 C d ~ — 7d (constant) 
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These forces act somewhere on the fin, and are signed again to match the vehicle frame, with 
7 > and a > 0. The summed forces on the body are thus: 

X = F a — \D\ cos a + \L\ sin a (116) 

- ~7a - Id + 7l« 2 

Z = F n + \L\ cos a + \D\ sin a 

— 7n« + 7l« + j D a. 

All of the forces are pushing the vehicle up. If we say that the fins act a distance If behind 
the temporary origin O, and that the moment carried by the fins themselves is very small 
(compared to the moment induced by Llf) the total moment is as follows: 

M = (- 7m + ln l n )a + ( 7£ + j D )l f a. (117) 
The moment about O vanishes if 



7m = iJn + IfilL + Id)- (118) 

The Munk moment j m opposes the aggregate effects of vorticity lift j n and the fins' lift and 
drag 7 £ + j D . With very large fins, this latter term is large, so that If might be very small; 
this is the case of AC moving aft toward the fins. A vehicle with excessively large fins will 
be difficult to turn and maneuver. 

Equation 118 contains two length measurements, referenced to an arbitrary body point O. 
To solve it explicitly, let lf n denote the (positive) distance that the fins are located behind 
F n ] this is likely a small number, since both effects usually act near the stern. We solve for 

If 

In + 7l + Id 

This is the distance that AC is located forward of the fins, and thus AC can be referenced 
to any other fixed point easily. Without fins, it can be recalled that 



i 1m 

In 

Hence, the fins act directly in the denominator to shorten If. Note that if the fins are located 
forward of the vortex shedding force F n , i.e., lf n < 0, If is reduced, but since AC is referenced 
to the fins, there is no net gain in stability. 

8.6 Aggregate Effects of Body and Fins 

Since all of the terms discussed so far have the same dependence on a, it is possible to group 
them into a condensed form. Setting F a and F n to account for the fuselage and fins, we have 
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X = F a ~- l - P C a A U 2 (120) 
Z = F n ~^pC' n A U 2 a 
M = -F n x AC ~-^C' n A U 2 x Ac a. 

8.7 Coefficients Z w , M w , Z q , and M q for a Slender Body 

The angle of attack a is related to the cross-body velocity w as follows: 

a = -tarT 1 ^) (121) 

W r tt 

~ — — lor U >> w. 
We can then write several linear hydrodynamic coefficients easily: 

Z w = -\pC' n A U (122) 

M w = l -pC' n A Ux Ac . 

The rotation of the vessel involves complex flow, which depends on both w and q, as well 
as their derivatives. To start, we consider the contribution of the fins only - slender body 
theory, introduced shortly, provides good results for the hull. The fin center of pressure 
is located a distance If aft of the body origin, and we assume that the vehicle is moving 
horizontally, with an instantaneous pitch angle of 9. The angle of attack seen by the fin is a 
combination of a part due to 9 and a part linear with q: 

a f ^-9+ l -g- (123) 
and so lateral force and moment derivatives (for the fin alone) emerge as 

Z q = - l -pC[A f Ul f (124) 
M q = - l -pC[A f Ul 2 f . 

9 SLENDER-BODY THEORY 

9.1 Introduction 



Consider a slender body with d « L, that is mostly straight. The body could be asymmetric 
in cross-section, or even flexible, but we require that the lateral variations are small and 
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smooth along the length. The idea of the slender-body theory, under these assumptions, 
is to think of the body as a longitudinal stack of thin sections, each having an easily- 
computed added mass. The effects are integrated along the length to approximate lift force 
and moment. Slender-body theory is accurate for small ratios d/L, except near the ends of 
the body. 

As one example, if the diameter of a body of revolution is d(s), then we can compute 5m a (x), 
where the nominal added mass value for a cylinder is 

Sm a = p^d 2 6x. (125) 

The added mass is equal to the mass of the water displaced by the cylinder. The equation 
above turns out to be a good approximation for a number of two-dimensional shapes, includ- 
ing flat plates and ellipses, if d is taken as the width dimension presented to the flow. Many 
formulas for added mass of two-dimensional sections, as well as for simple three-dimensional 
bodies, can be found in the books by Newman and Blevins. 



9.2 Kinematics Following the Fluid 

The added mass forces and moments derive from accelerations that fluid particles experience 
when they encounter the body. We use the notion of a fluid derivative for this purpose: the 
operator d/dt indicates a derivative taken in the frame of the passing particle, not the vehicle. 
Hence, this usage has an indirect connection with the derivative described in our previous 
discussion of rigid-body dynamics. 

For the purposes of explaining the theory, we will consider the two-dimensional heave/surge 
problem only. The local geometry is described by the location of the centerline; it has vertical 
location (in body coordinates) of Zb(x,t), and local angle a(x,t). The time- dependence 
indicates that the configuration is free to change with time, i.e., the body is flexible. Note 
that the curvilinear coordinate s is nearly equal to the body-reference (linear) coordinate x. 
The velocity of a fluid particle normal to the body at x is w n (t, x): 

w n = -^j- cos a — U sin a. (126) 

The first component is the time derivative in the body frame, and the second due to the 
deflection of the particle by the inclined body. If the body reference frame is rotated to the 
flow, that is, if w ^ 0, then dz^/dt will contain w. For small angles, sin a ~ tana = dzb/dx, 
and we can write 



dz b dz b 

The fluid derivative operator in action is as follows: 



dz b ( d d \ 



9.3 Derivative Following the Fluid 
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9.3 Derivative Following the Fluid 

A more formal derivation for the fluid derivative operator is quite simple. Let fi(x, t) represent 
some property of a fluid particle. 



^[li(x,t)] = Jim ^[fx(x + 5x,t + 8t) - fi(x,t)] 

_ u— 

dt dx 

The second equality can be verified using a Taylor series expansion of fi(x + Sx, t + St): 

u(x + 5x, t + St) = u(x, t) + — St + — Sx + h.o.t., 

ot ox 

and noting that Sx = —USt. The fluid is convected downstream with respect to the body. 
9.4 Differential Force on the Body 

If the local transverse velocity is w n (x, £), then the differential inertial force on the body here 
is the derivative (following the fluid) of the momentum: 

SF = -^[m a (x,t)w n (x,t)]Sx. (127) 

Note that we could here let the added mass vary with time also - this is the case of a changing 
cross- sect ion! The lateral velocity of the point Zb(x) in the body-reference frame is 

^ = w - xq, (128) 

such that 

w n = w — xq — Ua. (129) 

Taking the derivative, we have 



— = -(-- U- [m a (x,t)(w(t) - xq(t) - Ua(x,t))\ . 

We now restrict ourselves to a rigid body, so that neither m a nor a may change with time. 
S F c) 

— = m a (x)(-w + xq) + U — [m a (x)(w - xq - Ua)] . (130) 
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9.5 Total Force on a Vessel 

The net lift force on the body, computed with strip theory is 



/ SFdx (131) 

J XT 



where xt represents the coordinate of the tail, and xn is the coordinate of the nose. Ex- 
panding, we have 



rx N rx N Q 

Z = I m a (x) [— w + xq] dx + U / — — [m a (x)(w — xq — Ua)] dx 

J xt J xt OX 



= -m 33 w - m 35 q + Um a (x)(w - xq - Ua)\ x x= x x ™. 
We made use here of the added mass definitions 



rx N 

m 33 = / m a (x)dx 

J xt 

rx N 

m 35 — — I xm a {x)dx. 

J xt 



Additionally, for vessels with pointed noses and flat tails, the added mass m a at the nose is 
zero, so that a simpler form occurs: 

Z = —m 33 w — m 35 q — Um a (xT)(w — Xxq — Uo>(xt))- (132) 
In terms of the linear hydrodynamic derivatives, the strip theory thus provides 



Zy, = -m 33 

Zg = -m 35 

Z w = -Um a (x T ) 

Z q = Ux T m a (x T ) 

Z a ( XT ) = U 2 m a (x T ). 



It is interesting to note that both Z w and Z a ( XT ) depend on a nonzero base area. In general, 
however, potential flow estimates do not create lift (or drag) forces for a smooth body, so this 
should come as no surprise. The two terms are clearly related, since their difference depends 
only on how the body coordinate system is oriented to the flow. Another noteworthy fact is 
that the lift force depends only on a at the tail; a could take any value(s) along the body, 
with no effect on Z. 



9.6 Total Moment on a Vessel 
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9.6 Total Moment on a Vessel 

A similar procedure can be applied to the moment predictions from slender body theory 
(again for small a): 

fX N 

M — — xSFdx 

XN ( Q " 



rx N I Q (J \ 

= J x [-^- U -Q^j[ m a(x)(w-xq-Ua)]dx 

rxN rxN Q 

= I xm a (x)(w — xq)dx — U I x— [m a (x)(w — xq + Ua)] dx. 

J XT J XT OX 



fx N d 

'XT •> XT 

Then we make the further definition 



f XN 

m 55 — x m a (x)dx, 

J Xt 

(note that m 35 = m 53 ) and use integration by parts to obtain 



M = -m 35 w - m 55 q - Uxm a (x)(w - xq - Ua)\ x x= x x ^ + 

rXN 

U / m a (x)(w — xq — Ua)dx. 

J XT 



The integral above contains the product m a (x)a(x), which must be calculated if a changes 
along the length. For simplicity, we now assume that a is in fact constant on the length, 
leading to 

M = —rri35W — m,55q + UxTrri a (xT)(w — XTq — Ua) + 
Um 33 w + Um 35 q - U 2 m 33 a. 

Finally, the linear hydrodynamic moment derivatives are 



Myj = -m 35 

Mg = -m 55 

M w = Ux T m a (x T ) + Um 33 

M q = —Ux^m a {x T ) + Um 35 

M a = -U 2 x T m a (x T ) - U 2 m 33 . 

The derivative M w is closely-related to the Munk moment, whose linearization would provide 
M w = (m 33 — m u )U. The Munk moment (an exact result) may therefore be used to make 
a correction to the second term in the slender-body approximation above of M w . As with 
the lift force, M w and M a are closely related, depending only on the orientation of the body 
frame to the flow. 
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9.7 Relation to Wing Lift 

There is an important connection between the slender body theory terms involving added 
mass at the tail (m a (x T )), and low aspect-ratio wing theory. The lift force from the latter is 
of the form L = —\pUAC\w, where A = cs, the product of chord (long) and span (short). 
The lift coefficient slope is approximated by (Hoerner) 

C[ « \k{AR) : (133) 
where (AR) is the aspect ratio. Inserting this approximation into the lift formula, we obtain 

L = ~ps 2 Uw. (134) 

Now we look at a slender body approximation of the same force: The added mass at the tail 
is m a (xT) = ps 2 7r/4, and using the slender-body estimate for Z w , we calculate for lift: 

Z = —m a (x T )Uw 

= ——ps 2 Uw. 
4 

Slender-body theory is thus able to recover exactly the lift of a low-aspect ratio wing. Where 
does the slender-body predict the force will act? Recalling that M w = Um 33 + UxTm a (xT), 
and since 771,33 — for a front-back symmetric wing, the estimated lift force acts at the 
trailing edge. This location will tend to stabilize the wing, in the sense that it acts to orient 
the wing parallel to the incoming flow. 

9.8 Convention: Hydrodynamic Mass Matrix A 

Hydrodynamic derivatives that depend on accelerations are often written as components of 
a mass matrix A. By listing the body-referenced velocities in the order s — [u, v, w,p, q, r], 

— * 

we write (M + A)s = F, where M is the mass matrix of the material vessel and F is a 
generalized force. Therefore A S3 = —Z^, A 5t3 = — M^, and so on. 

10 PRACTICAL LIFT CALCULATIONS 

10.1 Characteristics of Lift-Producing Mechanisms 

At a small angle of attack, a slender body experiences transverse force due to: helical body 
vortices, the blunt trailing end, and fins. The helical body vortices are stable and symmetric 
in this condition, and are convected continuously into the wake. The low pressure associated 
with the vortices provides the suction force, usually toward the stern of the vehicle. The 
blunt trailing end induces lift as a product of added mass effects, and can be accurately 
modeled with slender body theory. A blunt trailing edge also induces some drag, which itself 
is stabilizing. The fins can often be properly modeled using experimental data parametrized 
with aspect ratio and several other geometric quantities. 



20.2 Jorgensen's Formulas 
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As the angle of attack becomes larger, the approximations in the fin and slender-body anal- 
ysis will break down. The helical vortices can become bigger while remaining stable, but 
eventually will split randomly. Some of it convects downstream, and the rest peels away 
from the body; this shedding is nonsymmetric, and greatly increases drag by widening the 
wake. In the limit of a 90° angle of attack, vorticity sheds as if from a bluff-body, and there 
is little axial convection. 



10.2 Jorgensen's Formulas 

There are some heuristic and theoretical formulas for predicting transverse force and moment 
on a body at various angles of attack, and we now present one of them due to Jorgensen. 
The formulas provide a good systematic procedure for design, and are best suited to vehicles 
with a blunt trailing edge. We call the area of the stern the base area. 
Let the body have length L, and reference area A r . This area could be the frontal projected 
area, the planform area, or the wetted area. The body travels at speed U, and angle of 
attack a. The normal force, axial force, and moment coefficients are defined as follows: 

Cn = T^nr (135) 



\pU 2 A r 

C A 



2' 

F A 



\p\J 2 A r 
" M x 



\pU 2 A r L' 

The moment M Xm is taken about a point x m , measured back from the nose; this location is 
arbitrary, and appears explicitly in the formula for C m . Jorgensen gives the coefficients as 
follows: 



Cn = -r- sin 2a cos ^ + -r-Cd n sin 2 a (136) 

C A = C Ao cos 2 a (137) 
V — Ab(L — x m ) . at A p ( x m — x^ 



C M = sin 2a cos - - C dn — y J sin a. (138) 

We have listed only the formulas for the special case of circular cross-section, although the 
complete equations do account for more complex shapes. Further, we have assumed that 
L » D, which is also not a constraint in the complete equations. The parameters used 
here are 

• A b : stern base area. A b = for a body that tapers to a point at the stern. 

• Cd n : crossflow drag coefficient; equivalent to that of an infinite circular cylinder. If 
"normal" Reynolds number Re n = U sin aD/u, 
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C d « 1.2, Re n < 3 x 10 5 

C d « 0.3, 3 x 10 5 < i?e n < 7 x 10 5 

C d « 0.6,7 x 10 5 < #e n . 



• ^4 P : planform area. 

• C J 4 o : axial drag at zero angle of attack, both frictional and form. C a — 0.002—0.006 for 
slender streamlined bodies, based on wetted surface area. It depends on Re = UL/v. 

• V: body volume. 

• x c : distance from the nose backwards to the center of the planform area. 

In the formula for normal force, we see that if ^4j, = 0, only drag forces act to create lift, 
through a sin 2 a-dependence. Similarly, the axial force is simply the zero-a result, modified 
by cos 2 a. In both cases, scaling of U 2 into the body principle directions is all that is required. 
There are several terms that match exactly the slender-body theory approximations for small 
a. These are the first term in the normal force (Cat), and the entire first term in the moment 
(Cm), whether or not Aj, = 0. Finally, we note that the second term in Cm disappears if 
x m = x c , i-e., if the moment is referenced to the center of the planform area. The idea here 
is that the fore and aft components of crossflow drag cancel out. 

The aerodynamic center (again referenced toward the stern, from the nose) can be found 
after the coefficients are computed: 



As written, the moment coefficient is negative if the moment destabilizes the body, while 
Cn is always positive. Thus, the moment seeks to move the AC forward on the body, but 
the effect is moderated by the lift force. 

10.3 Hoerner's Data: Notation 

An excellent reference for experimental data is the two-volume set by S. Hoerner. It con- 
tains a large amount of aerodynamic data from many different types of vehicles, wings, and 
other common engineering shapes. A few notations are used throughout the books, and are 
described here. 

First, dynamic pressure is given as q = \pU 2 , such that two typical body lift coefficients are: 



xac = x m + 




(139) 



C Y 



Y 




DLq 
Y 



D 2 q 



The first version uses the rectangular planform area as a reference, while the second uses the 
square frontal area. Hence, Cy = Cy d D/L. Two moment coefficients are: 



10.4 Slender-Body Theory vs. Experiment 
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r - N 



CjVi 



iVi 
LD 2 q 



where N is the moment taken about the body mid-point, and N\ is taken about the nose. 
Note that the reference area for moment is the square frontal projection, and the reference 
length is body length L. The following relation holds for these definitions 



c Nl - c N + -jC Yd . 

The lift and moment coefficients are strongly dependent on angle of attack; Hoerner uses 
the notation 

dC N 

db 

dC Nl 

db 

db ' 

and so on, where b is the angle of attack, usually in degrees. It follows from above that 

C n .b = C n b + Cydb/2- 



C n b — 
C n .b = 

Cyb = 



10.4 Slender-Body Theory vs. Experiment 

In an experiment, the net moment is measured, comprising both the destabilizing part due to 
the potential flow, and the stabilizing part due to vortex shedding or a blunt tail. Comparison 
of the measurements and the theory allows us to place the action point of the suction force. 
This section gives the formula for this location in Hoerner's notation, and gives two further 
examples of how well the slender-body theory matches experiments. 

For L/D > 6, the slender-body (pure added mass) estimates give C n b — —0.015/deg, acting 
to destabilize the vehicle. The value compares well with — 0.027/cfeg for a long cylinder 
and — 0.018 /deg for a long ellipsoid; it also reduces to —0.009 for L/D = 4. Note that the 
negative sign here is consistent with Hoerner's convention that destabilizing moments have 
negative sign. 

The experimental lift force is typically given by Cyb — 0.003 /deg; this acts at a point on 
the latter half of the vehicle, stabilizing the angle. Because this coefficient scales roughly 
with wetted area, proportional to LD, it changes little with L/D. It can be compared 
with a low-aspect ratio wing, which achieves an equivalent lift of n(AR)/2 = 0.0027 for 
(AR) = 10 ~ D/L. 

The point at which the viscous forces act can then be estimated as the following distance 
aft of the nose: 
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x 
L 



a 



n-b 

~C, 



a 



nb 



ydb 



(140) 



The calculation uses experimental values of Cydb an d C n .b, the moment slope referenced to 
the nose. In the table following are values from Hoerner (p. 13-2, Figure 2) for a symmetric 
and a blunt-ended body. 





symmetric 


blunt 


L/D 


6.7 


6.7 


C n b 


-0.012 


-0.012 


Cyb 


0.0031 


0.0037 


Cydb 


0.021 


0.025 


C n .b 


0.0012 (stable) 


0.0031 (stable) 


C n b 


-0.0093 (unstable) 


-0.0094 (unstable) 


x/L 


0.63 


0.60 



In comparing the two body shapes, we see that the moment at the nose is much more stable 
(positive) for the body with a blunt trailing edge. At the body midpoint, however, both 
vehicles are equally unstable. The blunt-tailed geometry has a much larger lift force, but it 
acts too close to the midpoint to add any stability there. 

The lift force dependence on the blunt tail is not difficult to see, using slender-body theory. 
Consider a body, with trailing edge radius r. The slender-body lift force associated with 
this end is simply the product of speed U and local added mass m a (xT) (in our previous 
notation). It comes out to be 



Z = -pU 2 (nr 2 )(2a), 



(141) 



such that the first term in parentheses is an effective area, and the second is a lift coefficient. 
With respect to the area itr 2 , the lift curve slope is therefore 2/rad. Expressed in terms of 
the Hoerner reference area D 2 , the equivalent lift coefficient is C y db = 0.0044/ deg, where we 
made the assumption here that 2r/D = 0.4 for the data in the table. This lift difference, 
due solely to the blunt end condition, is consistent with measurements. 



10.5 Slender-Body Approximation for Fin Lift 

Let us now consider two fins of span s each, acting at the tail end of the vehicle. This is 
the case if the vehicle body tapers to a point where the fins have their trailing edge. The 
slender-body approximation of lift as a result of blunt-end conditions is 

Z = ns 2 pU 2 a. (142) 

Letting the aspect ratio be (AR) = (2s) 2 /Af, where Af is the total area of the fin pair, 
substitution will give a lift curve slope of 



C[ = ^(AR). 
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This is known as Jones' formula, and is quite accurate for (AR) ~ 1. It is inadequate for 
higher-aspect ratio wings however, overestimating the lift by about 30% when (AR) = 2, 
and worsening further as (AR) grows. 



Vessels traveling at significant speed typically use rudders, elevators, and other streamlined 
control surfaces to maneuver. Their utility arises mainly from the high lift forces they can 
develop, with little drag penalty. Lift is always defined to act in a direction perpendicular 
to the flow, and drag in the same direction as the flow. 



A lifting surface is nominally an extrusion of a streamlined cross-section: the cross-section 
has a rounded leading edge, sharp trailing edge, and a smooth surface. The theory of lifting 
surfaces centers on the Kutta condition, which requires that fluid particle streamlines do not 
wrap around the trailing edge of the surface, but instead rejoin with streamlines from the 
other side of the wing at the trailing edge. This fact is true for a non-stalled surface at any 
angle of attack. 

Since the separation point on the front of the section rotates with the angle of attack, it is 
clear that the fluid must travel faster over one side of the surface than the other. The reduced 
Bernoulli pressure this induces can then be thought of as the lift-producing mechanism. More 
formally, lift arises from circulation T: 



and then L = —pUT. Circulation is the integral of velocity around the cross-section, and a 
lifting surface requires circulation in order to meet the Kutta condition. 

11.2 Three-Dimensional Effects: Finite Length 

Since all practical lifting surfaces have finite length, the flow near the ends may be three- 
dimensional. Prandtl's inviscid theory provides some insight. Since bound circulation cannot 
end abruptly at the wing end, it continues on in the fluid, leading to so-called wing-tip 
vortices. This continuation causes induced velocities at the tips, and some induced drag. 
Another description for the wing-tip vortices is that the pressure difference across the surface 
simply causes flow around the end. 

A critical parameter which governs the extent of three-dimensional effects is the aspect ratio: 



11 FINS AND LIFTING SURFACES 



11.1 Origin of Lift 




(143) 



AR 



span span 2 



(144) 



chord area 



The second representation is useful for non-rectangular control surfaces. The effective span 
is taken to be the length between the free ends of a symmetric wing. If the wing is attached 
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to a wall, the effective span is twice the physical value, by reflection, and in this case the 
effective aspect ratio is therefore twice the physical value. 

The aspect ratio is a strong determinant of wing performance: for a given angle of attack, a 
larger aspect ratio achieves a higher lift value, but also stalls earlier. 
Lift is written as 

L = \pU 2 AC h (145) 

where A is the single-side area of the surface. For angles of attack a below stall, the lift 
coefficient C t is nearly linear with a: C t = C[a, where C[ is called the lift coefficient slope, 
and has one empirical description 

C 'i = i | i 1 | i — (146) 

2na n(AR) 2vr(AR) 2 

where a is in radians, a ~ 0.90, and AR is the effective aspect ratio. When AR — > oo, the 
theoretical and maximum value for C[ is 2n. 

The lift generated on a surface is the result of a distributed pressure field; this fact creates 
both a net force and a net moment. A single equivalent force acts at a so-called center of 
action xa, which depends mainly on the aspect ratio. For high AR, xa — c/4, measured 
back from the front of the wing. For low AR, xa — c/2. 



11.3 Ring Fins 

Ring fins are useful when space allows, since they are omnidirectional, and structurally more 
robust than cantilevered plane surfaces. The effective aspect ratio for a ring of diameter d 
is given as 

Ad 

AR = — . (147) 

7TC 



The effective area of the ring is taken as 



A e = ^dc, (148) 



and we thus have L = \pU A e C[a, where one formula for C[ is 



C[ = . (149) 

U - 06 + n(AR) 
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12.1 Introduction 

We discuss in this section the nature of steady and unsteady propulsion. In many marine 
vessels and vehicles, an engine (diesel or gas turbine, say) or an electric motor drives the 



12.2 Steady Propulsion of Vessels 



51 



propeller through a linkage of shafts, reducers, and bearings, and the effects of each part 
are important in the response of the net system. Large, commercial surface vessels spend 
the vast majority of their time operating in open-water and at constant speed. In this case, 
steady propulsion conditions are generally optimized for fuel efficiency An approximation 
of the transient behavior of a system can be made using the quasi-static assumption. In the 
second section, we list several low-order models of thrusters, which have recently been used 
to model and simulate truly unsteady conditions. 

12.2 Steady Propulsion of Vessels 

The notation we will use is as given in Table 1, and there are two different flow conditions to 
consider. Self-propelled conditions refer to the propeller being installed and its propelling the 
vessel; there are no additional forces or moments on the vessel, such as would be caused by a 
towing bar or hawser. Furthermore, the flow around the hull interacts with the flow through 
the propeller. We use an sp subscript to indicate specifically self-propulsion conditions. 
Conversely, when the propeller is run in open water, i.e., not behind a hull, we use an o 
subscript; when the hull is towed with no propeller we use a t subscript. When subscripts 
are not used, generalization to either condition is implied. Finally, because of similitude 
(using diameter D in place of L when the propeller is involved), we do not distinguish 
between the magnitude of forces in model and full-scale vessels. 



Rsp 


N 


hull resistance under self-propulsion 


Rt 


N 


towed hull resistance (no propeller attached) 


T 


N 


thrust of the propeller 


n e 


Hz 


rotational speed of the engine 


Tim 


Hz 


maximum value of n e 


n p 


Hz 


rotational speed of the propeller 


A 




gear ratio 


Qe 


Nm 


engine torque 


Qp 


Nm 


propeller torque 


Vg 




gearbox efficiency 


Pe 


W 


engine power 


P P 


W 


propeller shaft power 


D 


m 


propeller diameter 


U 


m/s 


vessel speed 


u p 


m/s 


water speed seen at the propeller 


Qm 


Nm 


maximum engine torque 


f 


kg/s 


fuel rate (or energy rate in electric motor) 


fm 


kg/s 


maximum value of / 



Table 1: Nomenclature 
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12.2.1 Basic Characteristics 

In the steady state, force balance in self-propulsion requires that 

R SP = T sp . (150) 

The gear ratio A is usually large, indicating that the propeller turns much more slowly than 
the driving engine or motor. The following relations define the gearbox: 



n e = \n p (151) 

and power follows as P p = rj g P e , for any flow condition. We call J = U p /n p D the advance 
ratio of the prop when it is exposed to a water speed U p \ note that in the wake of the vessel, 
U p may not be the same as the speed of the vessel U. A propeller operating in open water 
can be characterized by two nondimensional parameters which are both functions of J: 



T 

K T = — ~~ ~ (thrust coefficient) (152) 

Kn = - (torque coefficient). (153) 

The open-water propeller efficiency can be written then as 

T oU ^ J(U)K T 

Vo 2nn p Q Po 2ttK q " U ° 4J 

This efficiency divides the useful thrust power by the shaft power. Thrust and torque co- 
efficients are typically nearly linear over a range of J, and therefore fit the approximate 
form: 



K T (J) = fc-foJ (155) 
k q( j ) = 7i - 72^- 

As written, the four coefficients [/3i, /3 2 , 7i, 72] are usually positive, as shown in the figure. 
We next introduce three factors useful for scaling and parameterizing our mathematical 
models: 

• U p — U(l — w); w is referred to as the wake fraction. A typical wake fraction of 0.1, 
for example, indicates that the incoming velocity seen by the propeller is only 90% of 
the vessel's speed. The propeller is operating in a wake. 

In practical terms, the wake fraction comes about this way: Suppose the open water 
thrust of a propeller is known at a given U and n p . Behind a vessel moving at speed 
U, and with the propeller spinning at the same n p , the prop creates some extra thrust. 
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Figure 4: Typical thrust and torque coefficients. 



w scales U at the prop and thus J; w is then chosen so that the open water thrust 
coefficient matches what is observed. The wake fraction can also be estimated by 
making direct velocity measurements behind the hull, with no propeller. 

• Rt = R sp (l — t). Often, a propeller will increase the resistance of the vessel by creating 
low-pressure on its intake side (near the hull), which makes R sp > R t . In this case, t 
is a small positive number, with 0.2 as a typical value, t is called the thrust deduction 
even though it is used to model resistance of the hull; it is obviously specific to both 
the hull and the propeller(s), and how they interact. 

The thrust deduction is particularly useful, and can be estimated from published values, 
if only the towed resistance of a hull is known. 

• Qp — VrQpsp- The rotative efficiency t]r, which may be greater than one, translates 
self-propelled torque to open water torque, for the same incident velocity U p , thrust 
T, and rotation rate n p . t]r is meant to account for spatial variations in the wake 
of the vessel which are not captured by the wake fraction, as well as the turbulence 
induced by the hull. Note that in comparison with the wake fraction, rotative efficiency 
equalizes torque instead of thrust. 

A common measure of efficiency, the quasi-propulsive efficiency, is based on the towed resis- 
tance, and the self-propelled torque. 



RtU 



(156) 



Vqp 



2nn p Q Psp 
T (l-t)U pVR 
2nn p (l - w)Q p 



Po 
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T a and Q Po are values for the inflow speed U p , and thus that rj Q is the open- water propeller 
efficiency at this speed. It follows that T (U P ) = T sp , which was used to complete the above 
equation. The quasi-propulsive efficiency can be greater than one, since it relies on the towed 
resistance and in general R t > R sp - The ratio (1 — t)/(l — w) is often called the hull efficiency, 
and we see that a small thrust deduction t and a large wake fraction w are beneficial effects, 
but which are in competition. A high rotative efficiency and open water propeller efficiency 
(at U p ) obviously contribute to an efficient overall system. 

12.2.2 Solution for Steady Conditions 

The linear form of K T and Kq (Equation 156) allows a closed-form solution for the steady- 
operating conditions. Suppose that the towed resistance is of the form 

Rt = \ P C r A w U\ (157) 

where C r is the resistance coefficient (which will generally depend on Re and Fr), and A w 
is the wetted area. Equating the self-propelled thrust and resistance then gives 



F sp Rsp 



Rt/il - 1) 

K T (J(U p ))pn 2 p D\l-t) = \ P C r A w U 2 



(Pi ~ (32J{U p ))pn 2 D A (l - t) = \ P C r A w S 

2 (1 — wy 



J(U P ) = . (158) 

The last equation predicts the steady-state advance ratio of the vessel, depending only on 
the propeller open characteristics, and on the hull. The vessel speed can be computed by 
recalling that J{U) = U/n p D and U p = U(l — w), but it is clear that we need now to find n p . 
This requires a torque equation, which necessitates a model of the drive engine or motor. 

12.2.3 Engine/Motor Models 

The torque-speed maps of many engines and motors fit the form 

Q e = QmF(f/f m ,n e /n m ), (159) 

where F() is the characteristic function. For example, gas turbines roughly fit the curves 
shown in the figure (Rubis). More specifically, if FQ has the form 
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Figure 5: Typical gas turbine engine torque-speed characteristic for increasing fuel rates 
fi, /2, h- 



F(f/f f m,n e /n m ) = - (a-f- + b) ^ + ( c-f- + d) (160) 



fm J f^m \ fn 

n m /X 



~ a i IT + a 2 



then a closed-form solution for n e (and thus n p ) can be found. The manipulations begin by 
equating the engine and propeller torque: 



Q Po (J(u P )) = vrQ P s P (J(u p )) 



pn 2 p D*K Q (J(U p )) = VRVgWe 
pnjD 5 (7i - ^2J{U P )) = r] R T]g\Q m F(f / 'f m ,n e /n m ) 



n l VRVg^Qm ( n p 

' ~ a ll TTT + a 2 



(n m /\y P D 5 ( 7l - 72 J(U p ))(n m /\y \ (n m /X) 



n p 



— ea\ + \J e 2 a\ + 4e«2 



{n m /X) 



(161) 



Note that the fuel rate enters through both ql\ and «2- 

The dynamic response of the coupled propulsion and ship systems, under the assumption of 
quasi-static propeller conditions, is given by 

(m + m a )u = T sp - R sp (162) 
lirlpiip = i~igXQ e Q PsP - 

Making the necessary substitutions creates a nonlinear model with / as the input; this is 
left as a problem for the reader. 
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12.3 Unsteady Propulsion Models 

When accurate positioning of the vehicle is critical, the quasi-static assumption used above 
does not suffice. Instead, the transient behavior of the propulsion system needs to be consid- 
ered. The problem of unsteady propulsion is still in development, although there have been 
some very successful models in recent years. It should be pointed out that the models de- 
scribed below all pertain to open-water conditions and electric motors, since the positioning 
problem has been central to bluff vehicles with multiple electric thrusters. 
We use the subscript m to denote a quantity in the motor, and p for the propeller. 

12.3.1 One-State Model: Yoerger et al. 

The torque equation at the propeller and the thrust relation are 

I p lo p = XQ m - K^ujp\u p \ (163) 
T = C t Up\up\. (164) 

where I p is the total (material plus fluid) inertia reflected to the prop;the propeller spins 
at uj p radians per second. The differential equation in uj p pits the torque delivered by the 
motor against a quadratic-drag type loss which depends on rotation speed. The thrust is 
then given as a static map directly from the rotation speed. 

This model requires the identification of three parameters: I p , K w , and Ct- It is a first-order, 
nonlinear, low-pass filter from Q m to T, whose bandwidth depends directly on Q m . 

12.3.2 Two-State Model: Healey et al. 

The two-state model includes the velocity of a mass of water moving in the vicinity of 
the blades. It can accommodate a tunnel around the propeller, which is very common in 
thrusters for positioning. The torque equation, similarly to the above, is referenced to the 
motor and given as 

I m u m = -K w u m + K V V - Qp/X. (165) 

Here, represents losses in the motor due to spinning (friction and resistive), and K v is 
the gain on the input voltage (so that the current amplifier is included in K v ). The second 
dynamic equation is for the fluid velocity at the propeller: 

pAL-filp = -pAA(3(U p -U)\U P -U\+T. (166) 

Here A is the disc area of the tunnel, or the propeller disk diameter if no tunnel exists. L 
denotes the length of the tunnel, and 7 is the effective added mass ratio. Together, pALj is 
the added mass that is accelerated by the blades; this mass is always nonzero, even if there 
is no tunnel. The parameter A/3 is called the differential momentum flux coefficient across 
the propeller; it may be on the order of 0.2 for propellers with tunnels, and up to 2.0 for 
open propellers. 
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The thrust and torque of the propeller are approximated using wing theory, which invokes 
lift and drag coefficients, as well as an effective angle of attack and the propeller pitch. 
However, these formulae are static maps, and therefore introduce no new dynamics. As 
with the one-state model of Yoerger et al., this version requires the identification of the 
various coefficients from experiments. This model has the advantage that it creates a thrust 
overshoot for a step input, which is in fact observed in experiments. 



Modern underwater vehicles and surface vessels are making increased use of electrical ac- 
tuators, for all range of tasks including weaponry, control surfaces, and main propulsion. 
This section gives a very brief introduction to the most prevalent electrical actuators: The 
DC motor, the AC induction motor, and the AC synchronous motor. For the latter two 
technologies, we consider the case of three-phase power, which is generally preferred over 
single-phase because of much higher power density; three-phase motors also have simpler 
starting conditions. AC motors are generally simpler in construction and more robust than 
DC motors, but at the cost of increased control complexity. 

This section provides working parametric models of these three motor types. As with gas 
turbines and diesel engines, the dynamic response of the actuator is quite fast compared to 
that of the system being controlled, say a submarine or surface vessel. Thus, we concentrate 
on portraying the quasi-static properties of the actuator - in particular, the torque/speed 
characteristics as a function of the control settings and electrical power applied. 
The discussion below on these various motors is generally invertible (at least for DC and AC 
synchronous devices) to cover both motors (electrical power in, mechanical power out) and 
generators (mechanical power in, electrical power out). We will only cover motors in this 
section, however; a thorough treatment of generators can be found in the references listed. 
The book by Bradley (19??) has been drawn on heavily in the following. 

13.1 Basic Relations 
13.1.1 Concepts 

First we need the notion of a magnetic flux, analagous to an electrical current, denoted $; 
a common unit is the Weber or Volt-second. The flux density 



is simply the flux per unit area, given in Teslas such that IT = lW/m 2 . Correponding to 
electrical field (Volts/m) is the magnetic field intensity H, in Amperes/meter: 



13 ELECTRIC MOTORS 



B = $/A 



(167) 



B $ 



(168) 



H = 
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/i 4/T x 10' 7 Henries /meter is the permeability of free space, and \x r is a (dimensionless) 
relative permeability. The product n \i r represents the real permeability of the material, 
and is thus the analog of electrical conductivity. A small area A or low relative permeability 
drives up the field intensity for a given flux 

13.1.2 Faraday's Law 

The voltage generated in a conductor experiencing a time rate of change in magnetic flux is 
given as 




(169) 



This voltage is commonly called the back-electromotive force or back-e.m.f., since it typically 
opposes the driving current; it is in fact a limiting factor in DC motors. 

13.1.3 Ampere's Law 

Current passing through a conductor in a closed loop generates a perpendicular magnetic 
field intensity given by 



An important point is that N circular wraps of the same conductor carrying current / induce 
the field H = ttDNI, where D is the diameter of the circle. 

13.1.4 Force 

Forces are generated from the orthogonal components of magnetic flux density B and current 
I: 



The units of this force is N/m, and so represents a distributed force on the conductor. 




(170) 



F 



I x B. 



(171) 



13.2 DC Motors 



The DC motor in its simplest form can be described by three relations: 



e a = K$uj 
V = e a + R a i a 
T = K$i a , 



where 



13.2 



DC Motors 
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• K is a constant of the motor 

• $ is the airgap magnetic flux per pole (Webers) 

• u) is the rotational speed of the motor (rad/s) 

• e a is the back-e.m.f. 

• V is the applied voltage 

• R a is the armature resistance on the rotor 

• i a is the current delivered to the armature on the rotor 

• T is developed torque 

The magnetic field in a typical motor is stationary (on the stator) and is created by perma- 
nent magnets or by coils, i.e., Faraday's law. Current is applied to the rotor armature through 
slip rings, and thus the force on each conductor in the armature is given by F = i a x B. 
Back-e.m.f. is created because the conductors in the rotor rotate through the stationary 
field, causing a relative rate of change of flux. The armature voltage loop contains the back- 
e.m.f. plus the resistive losses in the windings. As expected, torque scales with the product 
of magnetic flux and current. 

There are three main varieties of DC motors, all of which make use of the relations above. 
Speed control of the DC motor is primarily through the voltage V, either directly or through 
pulse- width modulation, but the stator flux could also be controlled in the shunt /independent 
configurations. 

13.2.1 Permanent Field Magnets 

Here, the magnetic field is created by permanent magnets arranged in the stator, imposing 
a steady $. The product K<& is generally written as k t , the torque constant of a DC motor, 
and has units of Nm/A. When SI units are used, k t also describes back-e.m.f.. The three 
basic relations are thus rewritten 



e a 
V 



T 



k t u 

&a F a i a 



which leads via substitution to 



1 



V 



RaT] 



h 



h J 



, or 



T 



■t 



[V 



k t uj] . 



■a 
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This result indicates that the torque developed scales linearly with the applied voltage, but 
that it also scales negatively with the motor speed. Hence, at the speed u> = V/K t , no torque 
is created. Additionally, the maximum torque is created at zero speed. 

Control of these motors is through the voltage V, or, more commonly, directly through 
current i a , which gives torque directly. 

13.2.2 Shunt or Independent Field Windings 

The field created by the stator can be strengthened by replacing the permanent magnets with 
electromagnets. The field windings are commonly placed in series with the rotor circuit, in 
parallel with it (shunt connection), or, they may be powered from a completely separate 
circuit. The latter two cases are effectively equivalent, in the sense that current and hence 
the field strength can be modulated easily, through a variable resistance in the shunt case. 
We have 



ou = 



v R " T 



with the important property that the second term in brackets is small due to the increased 
field strength, compared with the permanent magnet case above. Thus, the motor speed 
is effectively independent of torque, which makes these motor types ideal for self-regulation 
applications. At very high torques and currents, however, the total available flux will be 
reduced because of field armature reactance; the speed starts to degrade as shown in the 
figure. 

13.2.3 Series Windings 

When the field windings are arranged in series with the rotor circuit, the flux is 

$ = K s i a , 

where K s is a constant of the field winding. This additional connection requires 



V = e a + (R a + R s )i a ; 

the field winding brings a new resistance R s into the voltage loop. It follows through the 
substitutions that 



T = KK s i 2 a ^ 



I a = 



T 



v KK S 

V Ra + Rs 

to 



VKK7T KK S 



13.3 Three-Phase Synchronous Motor 



61 



The effects of resistance are usually quite small, so that the first term dominates, leading 
to a nonlinear torque/speed characteristic. The starting torque from this kind of motor is 
exceptionally high, and the series field winding finds wide application in railway locomotives. 
At the same time, it should be observed that under light loading, the series motor may well 
self-destruct since there is no intrinsic upper limit to speed! 




a) rob) roc) ro 

Figure 6: Torque-speed characteristics of three types of DC motors: a) permanent field 
magnets, b) shunt or independent field winding, c) series field winding. 



Some variations on the series and shunt connections are common, and referred to as com- 
pound DC motors. These achieve other torque/speed curves, including increasing torque 
with increasing speed, which can offset the speed droop due to field armature reaction ef- 
fects in the shunt motor. 

13.3 Three-Phase Synchronous Motor 

The rotor is either fitted with permanent magnets or supplied with DC current to create a 
static field on the rotor. The stator field windings are driven with three (balanced) phases 
of an AC supply, such that a moving field is created which rotates around the stator. The 
torque exerted on the rotor tries to align the two fields, and so the rotor follows the rotating 
stator field at the same speed. Note that if the rotor speed lags that of the stator field, there 
is no net torque; hence the name synchronous motor. 

A simple model of the synchronous motor is straightforward. As with the DC motors, the 
voltage loop equation for a single phase on the stator gives 

V = e a + ji a X, 

where V, e a , and i a are now phasors (magnitude of V and e a measured with respect to 
ground), j = \J— 1, and X is the reactance (armature and stator leakage) of the machine. 
Then, let (p denote the angle between i a and V. Equating the electrical (all three phases) 
and mechanical power gives 
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3Vi a cos (p = Tlo. 

Next, let 5 denote the angle between the phasors e a and V. It follows from the voltage loop 
equation that 

e a sin<5 — > 

p 3Ve a sinS 

, or 

2u X 

P Vabe a ,absm5 
2uj X 

where 

• p is the number of poles on the rotor; two poles means one north pole and one south 
pole, etc. 

• uj is both the rotational speed of the rotor, and the rotational speed of the stator field 

• Vab is the line-to-line voltage, equal to \/3V 

• e a ,ab is the line-to- line back-e.m.f., equal to V3e a 

• delta is the angle by which the rotor field lags the stator field (rad) 

• X is the synchronous reactance 

The torque scales with sin 8, and thus the rotor lags the stator field when the motor is 
powering; in a generator, the stator field lags the rotor. If the load torque exceeds the 
available torque, the synchronous motor can slip one or more poles, causing a large transient 
disturbance. 

Speed control of the three-phase synchronous machine is generally through the frequency of 
the three-phase power supply, uj, with the assumption that adequate voltage and current are 
available. 

13.4 Three-Phase Induction Motor 

Like the synchronous machine, the induction machine has windings on the stator to create 
a rotating magnetic field at frequency uj. Letting the rotor speed be uj r , we see immediately 
that if uj 7^ ui r , a potential field will be induced on any conductor on the rotor. In the case 
of a squirrel-cage rotor design, the rotor is made of conductor bars which are shorted out 
through rings at the ends, and hence the potential field will cause a real current flow. Torque 
is then generated through the familiar F = I x B relation. The fact of unequal field and rotor 
speeds in the induction motor is related to several unique effects, leading to torque-speed 
characteristic which differs significantly from both the DC the AC synchronous machines. 



i a cos 4>X = 
T = 

T = 



13.4 Three-Phase Induction Motor 
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First we define the slip ratio 



a slip ratio of zero means co = co r (and hence zero torque because the magnetic field seen by 
the rotor is constant) while a slip ratio of one implies the rotor is stopped. Most induction 
motors are designed to operate at a small positive slip ratio, say 0.1-0.2, for reasons described 
below. 

Next, since the magnetic flux lines pass through the rotor, the number of ampere-turns on 
stator and rotor is equivalent, that is, they form an ideal inductor: 



N r I r = NJ S . (173) 

We consider per-unit quantities from here on, for which we set N = N r = N s and hence 
I r = I s . If the stator flux at a particular location is $ s = $ sinu;t, the associated voltage is 
e s = NdQ /dt = N&°co cos cot. On the rotor, the same flux applies, but it rotates more slowly: 
$ r = cfrosinu;^ = $ D sin suit. Hence the rotor voltage is e r = Nd$ r /dt = NQ sco cos scot. 
Then the RMS voltage of the stator and rotor sides of the inductive coupling are related by 



The voltage seen at the stator scales inversely with the slip ratio, for a constant voltage at 
the rotor. In per-unit terms, the current in the rotor and stator are equivalent, and this then 
indicates that the rotor impedance, seen from the stator, also scales inversely with the slip 
ratio: 



Z rs = -[R r + jsX r ] 
s 

= — R r + jX r . 
s 

The factor of s in the rotor inductance occurs because the field seen by the rotor is actually 
rotating at so;. 

Next, we construct the (one phase) Thevenin equivalent circuit of the stator: it has a voltage 
source V t , and equivalent resistance R and inductance X. This is to be paired with the rotor 
resistance and inductance, reflected to the stator, giving the following current 



^J(R r /s + R) 2 + (X r + Xy' 

Finally, we need to express the torque/speed characteristic of the machine. The mechanical 
power is P m = (1 — s)coT, while the power delivered across the airgap is P gap = I 2 R r /s. The 
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actual power dissipated in the copper is related to the real rotor resistance: Pi oss = I 2 R r , 
and hence the mechanical power is P m = S(P gap — Pi oss ) = 3P gap (l — s). It follows that the 
efficiency of the motor is simply r\ — 1 — s. Combining the mechanical power with the torque 
equation gives 

P 3P 

rp 1 m _ 01 gap 

(1 — s)u> UJ 

3V t 2 R r 

suj[(R r /s + R) 2 + (X r + X) 2 }' 

Maximum torque is developed at a slight slippage, with decreased values at lower speeds. 
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Figure 7: Torque-speed characteristics of a typical three-phase induction motor. 
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Vehicles which are towed have some similarities to the vehicles that have been discussed so 
far. For example, towed vehicles are often streamlined, and usually need good directional 
stability. Some towed vehicles might have active lifting surfaces or thrusters for attitude 
control. On the other hand, if they are to be supported by a cable, towed vehicles may be 
quite heavy in water, and do not have to be self-propelled. The cable itself is an important 
factor in the behavior of the complete towed system, and in this section, we concentrate on 
cable mechanics more than vehicle characteristics, which can generally be handled with the 
same tools as other vehicles, i.e., slender-body theory, wing theory, linearization, etc.. Some 
basic guidelines for vehicle design are given at the end of this section. 

Modern cables can easily exceed 5000m in length, even a heavy steel cable with 2cm diameter. 
The cables are generally circular in cross section, and may carry power conductors and 
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multiple communication channels (fiber optic). The extreme L/D ratio for these cables 
obviates any bending stiffness effects. 

Cable systems come in a variety of configurations, and one main division may be made simply 
of the density of the cable. Light-tether systems are characterized by neutrally-buoyant (or 
nearly so) cables, with either a minimal vehicle at the end, as in a towed array, or a vehicle 
capable of maneuvering itself, such as a remotely-operated vehicle. The towed array is a 
relatively high-velocity system that nominally streams out horizontally behind the vessel. 
An ROV, on the other hand, operates at low speed, and must have large propulsors to 
control the tether if there are currents. Heavy systems, in contrast, employ a heavy cable 
and possibly a heavy weight; the rationale is that gravity will tend to keep the cable vertical 
and make the deployment robust against currents and towing speed. The heavy systems will 
generally transmit surface motions and tensions to the towed vehicle much more easily than 
light-tether systems. 

We will not discuss light systems specifically here, but rather look at heavy systems. Most 
of the analysis can be adapted to either case, however. 



14.1 Statics 
14.1.1 Force Balance 

For the purposes of deriving the static configuration of a cable in a flow, we assume for the 
moment that that it is inextensible. Tension and hydrostatic pressure will elongate a cable, 
but the effect is usually a small percentage of the total length. 

We employ the curvilinear axial coordinate s, which we take to be zero at the bottom end of 
the cable; upwards along the cable is the positive direction. The free-body diagram shown 
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has the following components: 

• W n : net in- water weight of the cable per unit length. 

• R n (s): external normal force, per unit length. 

• Rt(s): external tangential force, per unit length. 

• T(s): local tension. 

• 4>(s): local inclination angle. 

Force balance in the tangential and normal coordinates gives two coupled equations for T 
and (f): 

dT 

— = W n sm<P-R t (175) 

T^f- = W n cos (j) + R n . (176) 
ds 

The external forces are primarily fluid drag; the tangential drag is controlled by a frictional 
drag coefficient C t , and the normal drag scales with a crossflow drag coefficient C n . In both 
cases, the fluid velocity vector, U horizontal toward the left, is to be projected onto the 
relevant axes, leading to 

Rt = -\pC t dU 2 cos 2 <j> (177) 
Rn = ~ P C n dU 2 sin 2 (f>. (178) 

Note that we simplified the drag laws from a usual form v\v\ to v 2 , since as drawn, < < 
tt/2. 
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The equations for T and cf) can be integrated along the cable coordinate s to find the cable's 
static configuration. Two boundary conditions are needed, and the common case is that 
a force balance on the vehicle, dominated by drag, weight, and the cable tension, provides 
both T(0) and (f)(0). For example, a very heavy but low-drag vehicle will impose 0(0) ~ 7r/2, 
with T(0) equal to the in-water weight of the vehicle. 
With regard to Cartesian coordinates x, y, the cable configuration follows 



(IT 

^ = cos^ (179) 
ds 

^ = sin0. (180) 

as 

The simultaneous integration of all four equations (T, </>, x, y) defines the cable configuration, 
and current dependency may be included, say U is a function of y. 



14.1.2 Critical Angle 

For very deep systems, the total weight of cable will generally exceed that the vehicle. This 
gives rise to a configuration in which the cable is straight for a majority of its length, but 
turns as necessary at the vehicle end, to meet the bottom boundary condition. In the straight 
part of the cable, normal weight and drag components are equalized. The uniform angle is 
called the critical angle C , and can be approximated easily. Let the relative importance of 
weight be given as 
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so that the condition d(j)/ds — requires from the force balance 

S cos 4> c — — sin 2 C = 0. 

We are considering the case of < (fi c < n/2. Substituting sin 2 (fi c = 1 — cos 2 C , we solve a 
quadratic equation and keep only the positive solution: 



cos0 c = VS 2 + 1 - 1. (181) 

In the case of a very heavy cable, S is large, and the linear approximation of the square root 
Vl + e 1 + e/2 gives 

^ 1 _^ 

COS ffl c ~ — — >■ 

*■ - (182) 

For a very light cable, S is small; the same approximation gives 

cos 4> c ~ 1 — S — > (f) c ~ v^25 . 
The table below gives some results of the exact solution, and the approximations. 



5 


exact 


5 » 1 


5 « 1 


0.1 


0.44 




0.45 


0.2 


0.61 




0.63 


0.5 


0.91 


0.57 


1.00 


1.0 


1.14 


1.07 


1.41 


2.0 


1.33 


1.32 




5.0 


1.47 


1.47 





14.2 Linearized Dynamics 
14.2.1 Derivation 

The most direct procedure for deriving useful linear dynamic equations for a planar cable 
problem is to consider the total tension and angle as made up of static parts summed with 
dynamic parts: 



T(s,t) = T(s) + f(s,t) 
<l>(s,t) = $(s) + <j>(s,t). 
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We also write the axial deflection with respect to the static configuration as p(s, t), and the 
lateral deflection q(s, t). It follows that <j> = dq/ds. Now augment the two static configuration 
equations with inertial components: 

d 2 P df df Tir . .- ~. i ^ / ± d P \ 2 



(m + m a )^ = (f + f)^ + ^-iy n cos(0 + 0) 



. + 

2 



Here the material mass of the cable per unit length is m, and its transverse added mass 
is m a . Note that avoiding the drag law form v\v\ again, we have implicitly assumed that 
Ucos(f) > \dp/dt\ and Usm(f> > \dq/dt\. If it is not the case, say U = 0, then equivalent 
linearization can be used for the quadratic drag. 

Now we perform the trigonometry substitutions in the weight terms, let (j) ~ 4> for the 
calculation of drag, and substitute the constitutive (Hooke's) law 

df = EA &p 
ds ds 2 ' 

The static solution cancels out of both governing equations, and keeping only linear terms 
we obtain 

d 2 p d 2 p -dq -dp 

m— = ea ~q^ ~ W n cos (j)— - pC t dU cos (j)— 

, .d 2 q rf,d 2 q _. A dpdd> TTr . -dq 1TT . -dq 

(m + m a )^f = T^ + EA^^ + W n sm(f)^-pC n dUsm<t>-^. 

The axial dynamics (p) couples with the lateral equation through the weight term — W n cos (p(p. 
The lateral dynamics (q) couples with the axial through the term Td(p/ds. An additional 
weight term W n sin (fidq/ds also appears. The uncoupled dynamics are both in the form of 
damped wave equations 

d 2 p .dp ^ . d 2 p 

m W + bt m = EA dT 2 

, ,d 2 q . dq r^d 2( l ^ T ■ jdq 

(m + m a )— + b n - = T— + H/ n sm0-, 

where we made the substitution 6 t = pCtdU cos (p and b n = pC n dU sin 0. To a linear ap- 
proximation, the out-of-plane vibrations of a cable are also governed by the second equation 
above. 
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Because of light damping in the tangential direction, heavy cables easily transmit motions 
and tensions along their length, and can develop longitudinal resonant conditions (next 
section). In contrast, the lateral cable motions are heavily damped, such that disturbances 
only travel a few tens or hundreds of meters before they dissipate. The nature of the lateral 
response, in and out of the towing plane, is a very slow, damped nonlinear filter. High- 
frequency vessel motions in the horizontal plane are completely missed by the vehicle, while 
low-frequency motions occur sluggishly, and only after a significant delay time. 





axial 


lateral 


wave speed 


[ea 

V tn 


J T{S) 

y m+m a 




FAST 


SLOW 


natural frequency 


nir j EA 
L \ m 


oh\l T{L J 2) \ 

\ L y m+m a J 


mode shape 


sine/cosine 


Bessel function 


damping 


C t ^ 0(0.01) 


C n ~ 0(1) 


disturbances 






travel down 


YES 


NO 


cable? 







14.2.2 Damped Axial Motion 

Mode Shape. The axial direction is of particular interest, since it is lightly damped and 
forced by the heaving of vessels in seas. Consider a long cable governed by the damped wave 
equation 

mp + b t p = EAp" (183) 

We use over-dots to indicate time derivatives, and primes to indicate spatial derivatives. At 
the surface, we impose the motion 

p(L,t) = P cos cot, (184) 

while the towed vehicle, at the lower end, is an undamped mass responding to the local 
tension variations: 

EA ^A = M p(0 ,t). (185) 

These top and bottom behaviors comprise the boundary conditions for the wave equation. 
We let p(s, t) = p(s) coso;t, so that 
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This admits the solution p(s) = c\ cos ks + c 2 sin ks, where 



k = (187) 

Note that k is complex when b t ^ 0. The top and bottom boundary conditions give, respec- 
tively, 



P — c\ cos kL + c 2 sin kL 
= ci + £c 2 , 

where S = EAk/u 2 M. These can be combined to give the solution 

„ S cos ks — sin ks 
P=P 5^kL^kL - (188) 
In the case that M — > 0, the scalar S — > oo, simplifying the result to p = Pcos ks/ cos kL. 

Dynamic Tension. It is possible to compute the dynamic tension via T = EAp 1 . We 
obtain 

~ _ . „ , 5 sin ks + cos ks 

T =- EAPk ScoM - si , k L - < 189 > 
There are two dangerous situations: 

• The maximum tension is T + |T| and must be less than the working load of the cable. 
This is normally problematic at the top of the cable, where the static tension is highest. 

• If |T| > T, the cable will unload completely and then reload with extremely high 
impulsive forces. This is known as snap loading; it occurs primarily at the vehicle, 
where T is low. 



Natural Frequency. The natural frequency can be found by letting b t — 0, and investi- 
gating the singularity of p, for which 6 cos kL = sin kL. In general, kL « 1, but we find 



that a first-order approximation yields u = y EA/LM, which is only a correct answer if 
M » mL, i.e., the system is dominated by the vehicle mass. Some higher order terms need 
to be kept. We start with better approximations for sin() and cos(): 



2 J \ 6 J 

Employing the definition for 5, and recalling that u> 2 = k 2 EA/m, we arrive at 
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If we match up to second order in kL, then 



EA/L 
M + mL/2' 

This has the familiar form of the square root of a stiffness divided by a mass: the stiffness 
of the cable is EA/L, and the mass that is oscillating is M + mL/2. In very deep water, the 
effects of mL/2 dominate; if p c = m/A is the density of the cable, we have the approximation 





A few examples are given below for a steel cable with E = 200 x 10 9 Pa, and p c = 7000kg /m 3 . 
The natural frequencies near wave excitation at the surface vessel must be taken into account 
in any design or deployment. Even if a cable can withstand the effects of resonance, it may 
be undesirable to expose the vehicle to these motions. Some solutions in use today are: 
stable vessels (e.g., SWATH), heave compensation through an active crane, a clump weight 
below which a light cable is employed, and an S-shaped length of cable at the bottom formed 
with flotation balls. 



L = 500m 


u n = 15. Or ad /s 


1000m 


7.6rad/s 


2000m 


3.7rad/s 


5000m 


l.hrad/s 



14.3 Cable Strumming 

Cable strumming causes a host of problems, including obvious fatigue when the amplitudes 
and frequencies are high. The most noteworthy issue with towing is that the vibrations may 
cause the normal drag coefficient C n to increase dramatically - from about 1.2 for a non- 
oscillating cable to as high as 3.5. This drag penalty decreases the critical angle of towing, 
so that larger lengths of cable are needed to reach a given depth, and the towed system lags 
further and further behind the surface vessel. The static tension will increase accordingly as 
well. 

Strumming of cables is caused by the proximity of a preferred vortex formation frequency us 
to the natural frequency of the structure u> n . This latter frequency can be obtained as a zero 
of the lightly-damped Bessel function solution of the lateral dynamics equation above. The 
preferred frequency of vortex formation is given by the empirical relation u>s = 2nSU/d, 
where S is the Strouhal number, about 0.16-0.20 for a large range of Re. Strumming of 
amplitude d/2 or greater can occur for 0.6 < u)s/u) n < 2.0. The book by Blevins is a good 
general reference. 



14.4 Vehicle Design 
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14.4 Vehicle Design 

The physical layout of a towed vehicle is amenable to the analysis tools of self-propelled 
vehicles, with the main exceptions that the towpoint presents a large mean force as well 
as some disturbances, and that the vehicle can be quite heavy in water. Here are basic 
guidelines to be considered: 

1. The towpoint must be located above the vehicle center of in- water weight, for basic 
roll and pitch stability. 

2. The towpoint should be forward of the aerodynamic center, for towing stability reasons. 

3. The combined center of mass (material and added mass) should be longitudinally 
between the towpoint and the aerodynamic center, and nearer the towpoint. This will 
ensure that high-frequency disturbances do not induce excessive pitching. 

4. The towpoint should be longitudinally forward of the center of in-air weight, so that 
the vehicle enters the water fins first, and self-stabilizes with U > 0. 

5. The center of buoyancy should be behind the in-water center of weight, so that the 
vehicle pitches downward at small U, and hence the net lift force is downward, away 
from the surface. 

Meeting all of these criteria simultaneously is no small feat, and the performance of the 
device is very sensitive to small perturbations in the geometry. For this reason, full-scale 
experiments are commonly used in the design process. 



15 TRANSFER FUNCTIONS & STABILITY 

The reader is referred to Laplace Transforms in the section MATH FACTS for preliminary 
material on the Laplace transform. Partial fractions are presented here, in the context of 
control systems, as the fundamental link between pole locations and stability. 

15.1 Partial Fractions 

Solving linear time-invariant systems by the Laplace Transform method will generally create 
a signal containing the (factored) form 

Y( s ) = K(s + Zl )(s + z 2 )---(s + z m ) 
(s + Pl )(s + p 2 ) ■ ■ ■ (s + p n ) 

Although for the moment we are discussing the signal Y(s), later we will see that dynamic 
systems are described in the same format: in that case we call the impulse response G(s) 
a transfer function. A system transfer function is identical to its impulse response, since 
L(5(t)) = 1. 
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The constants —Zi are called the zeros of the transfer function or signal, and — pi are the 
poles. Viewed in the complex plane, it is clear that the magnitude of Y(s) will go to zero at 
the zeros, and to infinity at the poles. 

Partial fraction expansions alter the form of Y(s) so that the simple transform pairs can be 
used to find the time-domain output signals. We must have m < n; if this is not the case, 
then we have to divide the numerator by the denominator as necessary to find a simple form. 



15.2 Partial Fractions: Unique Poles 

Under the condition m < n, it is a fact that Y(s) is equivalent to 



Y(s) 



ai 



+ 



a 2 



+ 



(191) 



S+pi S+p 2 S +p n 

in the special case that all of the poles are unique and real. The coefficient a« is termed the 
residual associated with the z'th pole, and once all these are found it is a simple matter to 
go back to the transform table and look up the time-domain responses. 
How to find <2j? A simple rule applies: multiply the right-hand sides of the two equations 
above by (s +Pi), evaluate them at s — —pi, and solve for a iy the only one left. 



15.3 Example: Partial Fractions with Unique Real Poles 



G(s) 



s(s + 6) 



-2s 



(a + 4)(a-l) 

Since we have a pure delay and m = n, we can initially work with G(s) / se~ 2s . We have 



s + 6 ai a 2 

(s + 4)(s-l) = J+4 + 7^1' glVmg 



ai 

«2 



(s+6)(s+4) 
(«+4)(«-l) 

( s +6)(^-l) 



(»+4)(«-l) 



2 

5 



7 

5 



Thus 



L-\G{s)/se- 2s ) 



"5 6 



-it 



5 6 



g(t) = 5(t - 2) + |e- 4 (*- 2 ) + 

The impulse response is needed to account for the step change at t — 2. Note that in 
this example, we were able to apply the derivative operator s after expanding the partial 
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fractions. For cases where a second derivative must be taken, i.e., m > n + 1, special care 
should be used when accounting for the signal slope discontinuity at t — 0. The more 
traditional method, exemplified by Ogata, may prove easier to work through. 
The case of repeated real roots may be handled elegantly, but this condition rarely occurs 
in applications. 

15.4 Partial Fractions: Complex-Conjugate Poles 

A complex-conjugate pair of poles should be kept together, with the following procedure: 
employ the form 

Y(s)= f M + &2 - + -*»- + ..., (192) 
(s+ Pl )(s + p 2 ) s + p 3 

where Pl = p 2 (complex conjugate). As before, multiply through by (s + Pi)(s + p 2 ), and 
then evaluate at s = —p\. 

15.5 Example: Partial Fractions with Complex Poles 

G(s) - S + 1 - blS + &2 + — • 
s(s + j)(s-j) (s + j)(s-j) s 



S + V 



J s=-3 

i + i 



s + 1 



{s + j)(s-j) 



s=0 



= -hj + b 2 
= -1 



b 2 = 1; also 
= a 3 = 1. 



Working out the inverse transforms from the table of pairs, we have simply (noting that 

C = o) 



g(t) = — cost + sint + l(t). 



15.6 Stability in Linear Systems 

In linear systems, exponential stability occurs when all the real exponents of e are strictly 
negative. The signals decay within an exponential envelope. If one exponent is 0, the 
response never decays or grows in amplitude; this is called marginal stability. If at least one 
real exponent is positive, then one element of the response grows without bound, and the 
system is unstable. 
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15.7 Stability <^ Poles in LHP 



In the context of partial fraction expansions, the relationship between stability and pole 
locations is especially clear. The unit step function l(t) has a pole at zero, the exponential 
e -at a p G j e a ^ an( j go on ^jj Q £ other pairs exhibit the same property: A system 
is stable if and only if all of the poles occur in the left half of the complex plane. Marginally 
stable parts correlate with a zero real part, and unstable parts to a positive real part. 

15.8 General Stability 

There are two definitions, which apply to systems with input u(t) and output y(t). 

1. Exponential. If u(t) = and 7/(0) = y Q , then \yi(t)\ < ae" 7 *, for finite a and 7 > 0. 
The output asymptotically approaches zero, within a decaying exponential envelope. 

2. Bounded-Input Bounded-Output (BIBO). If 7/(0) = 0, and \ fi(t)\ < 7,7 > and 
finite, then \yi(t)\ < a, a > and finite. 

In linear time-invariant systems, the two definitions are identical. Exponential stability is 
easy to check for linear systems, but for nonlinear systems, BIBO stability is usually easier 
to achieve. 



16.1 Introduction 

16.1.1 Plants, Inputs, and Outputs 

Controller design is about creating dynamic systems that behave in useful ways. Many 
target systems are physical; we employ controllers to steer ships, fly jets, position electric 
motors and hydraulic actuators, and distill alcohol. Controllers are also applied in macro- 
economics and many other important, non-physical systems. It is the fundamental concept 
of controller design that a set of input variables acts through a given "plant" to create an 
output. Feedback control then uses sensed plant outputs to apply corrective inputs: 
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Plant 



Inputs 



Outputs 



Sensors 



Jet aircraft 
Marine vessel 
Hydraulic robot 
U.S. economy 
Nuclear reactor 



fed interest rate, etc. 
cooling, neutron flux 



elevator, rudder, etc. 
rudder angle 



valve position 



tip position 
prosperity 
power level 



altitude, hdg 
heading 



altimeter, GPS 
gyrocompass 
joint angle 
inflation, Ml 
temp., pressure 



16.2 Representing Linear Systems 
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16.1.2 The Need for Modeling 

Effective control system design usually benefits from an accurate model of the plant, although 
it must be noted that many industrial controllers can be tuned up satisfactorily with no 
knowledge of the plant. Ziegler and Nichols, for example, developed a general recipe which 
we detail later. In any event, plant models simply do not match real- world systems exactly; 
we can only hope to capture the basic components in the form of differential or integro- 
differential equations. 

Beyond prediction of plant behavior based on physics, the process of system identification 
generates a plant model from data. The process is often problematic, however, since the 
measured response could be corrupted by sensor noise or physical disturbances in the system 
which cause it to behave in unpredictable ways. At some frequency high enough, most 
systems exhibit effects that are difficult to model or reproduce, and this is a limit to controller 
performance. 

16.1.3 Nonlinear Control 

The bulk of this subject is taught using the tools of linear systems analysis. The main 
reason for this restriction is that nonlinear systems are difficult to model, difficult to design 
controllers for, and difficult overall! Within the paradigm of linear systems, there are many 
sets of powerful tools available. The reader interested in nonlinear control is referred to the 
book by Slotine and Li (1991). 

16.2 Representing Linear Systems 

Except for the most heuristic methods of tuning up simple systems, control system design 
depends on a model of the plant. The transfer function description of linear systems has 
already been described in the discussion of the Laplace transform. The state-space form is 
an entirely equivalent time-domain representation that makes a clean extension to systems 
with multiple inputs and multiple outputs, and opens the way to standard tools from linear 
algebra. 

16.2.1 Standard State-Space Form 

We write a linear system in a state-space form as follows 

x = Ax + Bu + Gw (193) 
y = Cx + Du + v 

where 

• x is a state vector, with as many elements as there are orders in the governing differ- 
ential equations. 

• A is a matrix mapping x to its derivative; A captures the natural dynamics of the 
system without external inputs. 
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• B is an input gain matrix for the control input u. 

• G is a gain matrix for unknown disturbance w; w drives the state just like the control 
u. 

• y is the observation vector, comprised mainly of a linear combination of states Cx 
(where C is a matrix). 

• Du is a direct map from input to output (usually zero for physical systems). 

• v is an unknown sensor noise which corrupts the measurement. 




16.2.2 Converting a State-Space Model into a Transfer Function 

There are a number of canonical state-space forms available, which can create the same 
transfer function. In the case of no disturbances or noise, the transfer function (or transfer 
matrix) can be written as 

G(s) = 44 = Chl-A^B + D, (194) 
u(s) 

where / is the identity matrix with the same size as A. A similar equation holds for y(s)/w(s), 
and clearly y(s)/v(s) = I. 

16.2.3 Converting a Transfer Function into a State-Space Model 

It may be possible to write the corresponding differential equation along one row of the state 
vector, and then cascade derivatives. For example, consider the following system: 



my"(t) + by' {t) + ky(t) = u'(t) + u(t) (mass-spring-dashpot) 

v ' ms 2 + bs + k 
Setting x = [y 1 , y] T , we obtain the system 
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dx 
~dt 

y 



—b/m 
1 

[1 l]x 



-k/m 




x + 



1/m 




u 



Note specifically that dx-ijdt = X\, leading to an entry of 1 in the off-diagonal of the second 
row in A. Entries in the C-matrix are easy to write in this case because of linearity; the 
system response to v! is the same as the derivative of the system response to u. 

16.3 PID Controllers 

The most common type of industrial controller is the proportional-integral-derivative (PID) 
design. If u is the output from the controller, and e is the error signal it receives, this control 
law has the form 



coo 



u(t) 

E{8) 



k p e(t) + ki j e(r)dr + k d e'(t), 

J 



kn 



k d s 



l + — +T d S 
TiS 



(195) 



where the last line is written using the conventions of one overall gain k p , plus a time 
characteristic to the integral part (r«) and and time characteristic to the derivative part (r^). 
In words, the proportional part of this control law will create a control action that scales 
linearly with the error - we often think of this as a spring-like action. The integrator 
is accumulating the error signal over time, and so the control action from this part will 
continue to grow as long as an error exists. Finally, the derivative action scales with the 
derivative of the error. The controller will retard motion toward zero error, which helps to 
reduce overshoot. 

The common variations are: P, PD, PI, PID. 
16.4 Example: PID Control 

Consider the case of a mass (m) sliding on a frictionless table. It has a perfect thruster 
that generates force u(t), but is also subject to an unknown disturbance d(t). If the linear 
position of the mass is y(t), and it is perfectly measured, we have the plant 



my" {t) = u(t) + d(t). 

Suppose that the desired condition is simply y(t) = 0, with initial conditions y(0) = y Q and 
2/(0) =0. 
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16.4.1 Proportional Only 

A proportional controller alone invokes the control law u(t) = —k p y(t), so that the closed- 
loop dynamics follow 



my" (t) = -k p y(t) + d(t). 

In the absence of d(t), we see that y(t) = y cos \f^t, a marginally stable response that is 
undesirable. 

16.4.2 Proportional-Derivative Only 

Let u(t) = —k p y(t) — k d y'(t), and it follows that 



my"(t) = -k p y{t) - k d y'(t) + d(t). 

The system now resembles a second-order mass-spring-dashpot system where k p plays the 
part of the spring, and k d the part of the dashpot. With an excessively large value for 
k d , the system would be overdamped and very slow to respond to any command. In most 
applications, a small amount of overshoot is employed because the response time is shorter. 
The k d value for critical damping in this example is 2^Jmk p , and so the rule is k d < 2^Jmk p . 
The result, easily found using the Laplace transform, is 



y(t) = y e *» 



COS (jJ d t + 



2mu d 



sin u) d t 



where cu d = yAmkp — k\j2m. This response is exponentially stable as desired. Note that if 
the mass had a very large amount of natural damping, a negative k d could be used to cancel 
some of its effect and speed up the system response. 

Now consider what happens if d(t) has a constant bias d Q : it balances exactly the proportional 
control part, eventually settling out at y(t = oo) = d /k p . To achieve good rejection of d 
with a PD controller, we would need to set k p very large. However, very large values of k p 
will also drive the resonant frequency u> d up, which is unacceptable. 



16.4.3 Proportional-Integral-Derivative 

Now let u(t) = —k p y{t) — ki J^y^dr — k d y'(t): we have 



my" (t) = -k p y(t) - h Cy^dr - k d y'(t) + d{t). 

Jo 

The control system has now created a third-order closed-loop response. If d(t) = d , a time 
derivative leads to 



my"'(t) + k p y'{t) + k iV (t) + k d y"(t) = 0, 
so that y(t = oo) = 0, as desired, provided the roots are stable. 



16.5 Heuristic Tuning 
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16.5 Heuristic Tuning 

For many practical systems, tuning of a PID controller may proceed without any system 
model. This is especially pertinent for plants which are open-loop stable, and can be safely 
tested with varying controllers. One useful approach is due to Ziegler and Nichols (e.g., 
Belanger,1995), which transforms the basic characteristics of a step response (e.g., the input 
is 1(£)) into a reasonable PID design. The idea is to approximate the response curve by a 
first-order lag (gain k and time constant r) and a pure delay T: 

hp-Ts 

G(s) ~ — (196) 

w rs + 1 v ' 

The following rules apply only if the plant contains no dominating, lightly-damped complex 
poles, and has no poles at the origin: 



p 


k p 


= l.Or/T 






PI 


k p 


= 0.9r/T 


ki = 0.27r/T 2 
h = O.QOt/T 2 




PID 


k p 


= 1.2r/T 


k d = 0.60r 



Note that if no pure time delay exists (T = 0), this recipe suggests the proportional gain can 
become arbitrarily high! Any characteristic other than a true first-order lag would therefore 
be expected to cause a measurable delay. 

16.6 Block Diagrams of Systems 
16.6.1 Fundamental Feedback Loop 

The topology of a feedback system can be represented graphically by considering each dy- 
namical system element to reside within a box, having an input line and an output line. For 
example, the plant used above (a simple mass) has transfer function P(s) = 1/ms 2 , which 
relates the input, force u(s), into the output, position y(s). In turn, the PD-controller has 
transfer function C(s) = k p + k d s; its input is the error signal E(s) = —y(s), and its output 
is force u(s). The feedback loop in block diagram form is shown below. 





as) 


u(s) 


P(s) 


y(s) 




^- 







16.6.2 Block Diagrams: General Case 

The simple feedback system above is augmented in practice by three external inputs. The 
first is a process disturbance, which can be taken to act at the input of the physical plant, 
or at the output. In the former case, it is additive with the control action, and so has 
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some physical meaning. In the second case, the disturbance has the same units as the plant 
output. 

Another external input is the reference command or setpoint, used to create a more general 
error signal e(s) = r(s) — y(s). Note that the feedback loop, in trying to force e(s) to zero, 
will necessarily make y(s) approximate r(s). 

The final input is sensor noise, which usually corrupts the feedback signal y(s), causing 
some error in the evaluation of e(s), and so on. Sensors with very poor noise properties can 
ruin the performance of a control system, no matter how perfectly understood are the other 
components. 



+ ^ e 









c 




p 





+ 



16.6.3 Primary Transfer Functions 

Some algebra shows that 



e 


1 


r 


1 + PC 


y 


PC 


r 


1 + PC 


u 


c 
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e/r = S relates the reference input and noise to the error, and is known as the sensitivity 
function. We would generally like S to be small at low frequencies, so that the tracking 
error there is small, y/r = T is called the complementary sensitivity function. Note that 
S + T = 1, implying that these two functions must always trade off; they cannot both be 
small or large at the same time. Other systems we encounter again later are the (forward) 
loop transfer function PC, the loop transfer function broken between C and P: CP, and 



e_ -P 

~d u ~ 1 + PC 

y_ P 

d u ' 1 + PC 

u_ -CP 

~d u ~ 1 + CP 

e -1 
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y_ i 

d y 1 + PC 

v_ -C 

dy~ ~ 1 + CP 

e -1 

n ~ 1 + PC 
y -PC 

n 1 + PC 
u -C 
n ~ 1 + CP 

If the disturbance is taken at the plant output, then the three functions S, T, and U (control 
action) completely describe the system. This will in fact be the procedure when we address 
loopshaping. 

17 MODAL ANALYSIS 

17.1 Introduction 

The evolution of states in a linear system occurs through independent modes, which can be 
driven by external inputs, and observed through plant output. This section provides the 
basis for modal analysis of systems. Throughout, we use the state-space description of a 
system with D = 0: 

x = Ax + Bu 
y = Cx. 

17.2 Matrix Exponential 
17.2.1 Definition 

In the instance of an unforced response to initial conditions, consider the system 

x = Ax, x{t = 0) = x- 

In the scalar case, the response is x(t) = xe a *, giving a decaying exponential if a < 0. The 
same notation holds for the case of a vector x, and matrix A: 

x(t) = e At x, where 
e At = I + At+ { -^ + ... 

e At is usually called the matrix exponential. 



= S 

= -u 
= -s 

= -T 
= -U. 
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17.2.2 Modal Canonical Form 

Introductory material on the eigenvalue problem and modal decomposition can be found in 
the MATH FACTS section. This modal decomposition of A leads to a very useful state-space 
representation. Namely, since A = VAV" 1 , a transformation of state variables can be made, 
x = Vz, leading to 



'z = Az + V^Bu (197) 
y = CVz. 

This is called the modal canonical form, since the states are simply the modal amplitudes. 
These states are uncoupled in A, but may be coupled through the input (V^B) and output 
(CV) mappings. The modal form is numerically robust for computations. 

17.2.3 Modal Decomposition of Response 

Now we are ready to look at the matrix exponential e At in terms of its constituent modes. 
Employing the above form for A, we find that 



e At = I + At+ { -^ + ... 
= Ve At W T 

n 
i=l 

In terms of the response to an initial condition x, we have 



i=i 

The product w[ x * s a scalar, the projection of the initial conditions onto the i'th mode. 
If x is perpendicular to Wi, then the product is zero and the i'th mode does not respond. 
Otherwise, the i'th mode does participate in the response. The projection of the i'th mode 
onto the states x is through the right eigenvector t^. 

For stability of the system, the eigenvalues of A, that is, Aj, must have negative real parts; 
they are in fact the poles of the equivalent transfer function description. 

17.3 Forced Response and Controllability 

Now consider the system with an external input u: 



17.4 Plant Output and Observability 
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x = Ax + Bu, x(t = 0) = x- 

Taking the Laplace transform of the system, taking into account the initial condition for the 
derivative, we have 

sx(s) — x — Ax(s) + Bu(s) — > 

x(s) = (sI-A^x+isI-Ay'Buis). 

Thus (si — A)^ 1 can be recognized as the Laplace transform of the matrix exponential e . 
In the time domain, the second term then has the form of a convolution of the matrix 
exponential and the net input Bu: 

x(t) = t e A{t - T) Bu(T)dr 
Jo 

= E t e Xi{t ~ T) ViwjBu(T)dT. 
i=i Jo 

Suppose now that there are m inputs, such that B — [&i, 62, • • • , b m ]. Then some rearrange- 
ment will give 

i=i k=i Jo 

The product wfbk, a scalar, represents the projection of the fc'th control channel onto the 
i'th mode. We say that the i'th mode is controllable from the fc'th input if the product 
is nonzero. If a given mode % has wj bk = for all input channels k, then the mode is 
uncontrollable. 

In normal applications, controllability for the entire system is checked using the following 
test: Construct the so-called controllability matrix: 

M c = [B, AB, • • • , A n ~ l B}. (198) 

This matrix has size n x (nm), where m is the number of input channels. If M c has rank n, 
then the system is controllable, i.e., all modes are controllable. 

17.4 Plant Output and Observability 

We now turn to a related question: can the complete state vector of the system be observed 
given only the output measurements y, and the known control ul The response due to 
the external input is easy to compute deterministically, through the convolution integral. 
Consider the part due to initial conditions x- We found above 
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n 



m = E 



i=l 



The observation is y = Cx (r channels of output), and writing 



C = 




the fc'th channel of the output is 



n 



i=l 



The i'th mode is observable in the fc'th output if the product c^Vi ^ 0. We say that a 
system is observable if every mode can be seen in at least one output channel. The usual 
test for system observability requires computation of the observability matrix: 



18.1 Introduction 

This section formalizes the notion of loopshaping for linear control system design. The 
loopshaping approach is inherently two-fold. First, we shape the open- loop transfer function 
(or matrix) P(s)C(s), to meet performance and robustness specifications. Once this is done, 
then the compensator must be computed, from from knowing the nominal product P(s)C(s), 
and the nominal plant P(s). 

Most of the analysis here is given for single- input, single-output systems, but the link to 
multivariable control is not too difficult. In particular, absolute values of transfer functions 
are replaced with the maximum singular values of transfer matrices. Design based on singular 
values is the idea of L 2 -control, or LQG/LTR, to be presented in the next lectures. 

18.2 Roots of Stability — Nyquist Criterion 

We consider the SISO feedback system with reference trajectory r(s) and plant output y(s), 
as given previously. The tracking error signal is defined as e(s) = r(s) — y(s), thus forming 
the negative feedback loop. The sensitivity function is written as 




(199) 
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18.2 Roots of Stability - Nyquist Criterion 
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S(s) 



e(s) _ 1 



r(s) l + P(s)C(s)' 



where P(s) represents the plant transfer function, and C(s) the compensator. The closed- 
loop characteristic equation, whose roots are the poles of the closed-loop system, is 1 + 
P(s)C(s) = 0, equivalent to P(s)C_(s) + P(s)C(s) = 0, where the underline and overline 
denote the denominator and numerator, respectively. The Nyquist criterion allows us to 
assess the stability properties of a system based on P(s)C(s) only. This method for design 
involves plotting the complex loci of P(s)C(s) for the range u> = [—00,00]. There is no 
explicit calculation of the closed-loop poles, and in this sense the design approach is quite 
different from the root-locus method (see Ogata). 

18.2.1 Mapping Theorem 

We impose a reasonable assumption from the outset: The number of poles in P(s)C(s) 
exceeds the number of zeros. It is a reasonable constraint because otherwise the loop transfer 
function could pass signals with infinitely high frequency. In the case of a PID controller 
(two zeros) and a second-order zero- less plant, this constraint can be easily met by adding 
a high-frequency rolloff to the compensator, the equivalent of low-pass filtering the error 
signal. 

Let F(s) = 1 + P(s)C(s). The heart of the Nyquist analysis is the mapping theorem, which 
answers the following question: How do paths in the s-plane map into paths in the F-plane? 
We limit ourselves to closed, clockwise(CW) paths in the s-plane, and the remarkable result 
of the mapping theorem is 

Every zero of F(s) enclosed in the s-plane generates exactly one CW encirclement of the 
origin in the F(s) -plane. Conversely, every pole of F(s) enclosed in the s-plane generates 
exactly one CCW encirclement of the origin in the F(s) -plane. Since CW and CCW encir- 
clements of the origin may cancel, the relation is often written Z — P = CW. 
The trick now is to make the trajectory in the s-plane enclose all unstable poles, i.e., the 
path encloses the entire right-half plane, moving up the imaginary axis, and then proceeding 
to the right at an arbitrarily large radius, back to the negative imaginary axis. 
Since the zeros of F(s) are in fact the poles of the closed-loop transfer function, e.g., S(s), 
stability requires that there are no zeros of F(s) in the right-half s-plane. This leads to a 
slightly shorter form of the above relation: 



In words, stability requires that the number of unstable poles in F(s) is equal to the number 
of CCW encirclements of the origin, as s sweeps around the entire right-half s-plane. 

18.2.2 Nyquist Criterion 

The Nyquist criterion now follows from one translation. Namely, encirclements of the origin 
by F(s) are equivalent to encirclements of the point (—1 + Oj) by F(s) — 1, or P(s)C(s). 



p = CCW. 



(200) 
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Then the stability criterion can be cast in terms of the unstable poles of P(s)C(s), instead 
of those of F(s): 

P = CCW i — > closed-loop stability (201) 

This is in fact the complete Nyquist criterion for stability. It is a necessary and sufficient 
condition that the number of unstable poles in the loop transfer function P(s)C(s) must be 
matched by an equal number of CCW encirclements of the critical point (—1 + Oj). 
There are several details to keep in mind when making Nyquist plots: 

• If neither the plant nor the controller have unstable modes, then the loci of P(s)C(s) 
must not encircle the critical point at all. 

• Because the path taken in the s-plane includes negative frequencies (i.e., the nega- 
tive imaginary axis), the loci of P(s)C(s) occur as complex conjugates - the plot is 
symmetric about the real axis. 

• The requirement that the number of poles in P(s)C(s) exceeds the number of zeros 
means that at high frequencies, P(s)C(s) always decays such that the loci go to the 
origin. 

• For the multivariable (MIMO) case, the procedure of looking at individual Nyquist 
plots for each element of a transfer matrix is unreliable and outdated. Referring to 
the multivariable definition of S(s), we should count the encirclements for the function 
[det(I + P(s)C(s)) — 1] instead of P(s)C(s). The use of gain and phase margin in 
design is similar to the SISO case. 

18.2.3 Robustness on the Nyquist Plot 

The question of robustness in the presence of modelling errors is central to control system 
design. There are two natural measures of robustness for the Nyquist plot, each having a 
very clear graphical representation. The loci need to stay away from the critical point; how 
close the loci come to it can be expressed in terms of magnitude and angle. 

• When the angle of P(s)C(s) is —180°, the magnitude \P(s)C(s)\ should not be near 
one. 

• When the magnitude \P(s)C(s) \ = 1, its angle should not be —180°. 

These notions lead to definition of the gain margin k 9 and phase margin 7 for a design. As 
the figure shows, the definition of k g is different for stable and unstable P(s)C(s). Rules of 
thumb are as follows. For a stable plant, k g > 2 and 7 > 30°; for an unstable plant, k g < 0.5 
and 7 > 30°. 



18.3 Design for Nominal Performance 
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18.3 Design for Nominal Performance 

Performance requirements of a feedback controller, using the nominal plant model, can be 
cast in terms of the Nyquist plot. We restrict the discussion to the scalar case in the following 
sections. 

Since the sensitivity function maps reference input r(s) to tracking error e(s), we know 
that (S^s)! should be small at low frequencies. For example, if one-percent tracking is to 
be maintained for all frequencies below uj = A, then \S(s)\ < 0.01, \/u < X. This can be 
formalized by writing 



where Wi(s) is a stable weighting function of frequency. To force S(s) to be small at low 
u>, Wi{s) should be large in the same range. The requirement |H / i(s)S'(s)| < 1 is equivalent 
to |Wi(s)| < |1 + P(s)C(s)\, and this latter condition can be interpreted as: The loci of 
P(s)C(s) must stay outside the disk of radius Wi(s), which is to be centered on the critical 
point (— 1+0 j). The disk is to be quite large, possibly infinitely large, at the lower frequencies. 

18.4 Design for Robustness 

It is a ubiquitous observation that models of plants degrade with increasing frequency. For 
example, the DC gain and slow, lightly-damped modes or zeros are easy to observe, but 
higher-frequency components in the response may be hard to capture or even excite repeat- 
ably. Higher- frequency behavior may have more nonlinear properties as well. 
The effects of modeling uncertainty can be considered to enter the nominal feedback system 
as a disturbance at the plant output, d y . One of the most useful descriptions of model 
uncertainty is the multiplicative uncertainty: 




(202) 



P(s) = (l + A(s)W 2 (s))P(s). 



(203) 
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Here, P(s) represents the nominal plant model used in the design of the control loop, 
and P(s) is the actual, perturbed plant. The perturbation is of the multiplicative type, 
A(s)W2{s)P(s), where A(s) is an unknown but stable function of frequency for which 
| A(s)| < 1. The weighting function W 2 (s) scales A(s) with frequency; W 2 (s) should be 
growing with increasing frequency, since the uncertainty grows. However, W 2 (s) should not 
grow any faster than necessary, since it will turn out to be at the cost of nominal performance. 
In the scalar case, the weight can be estimated as follows: since P/P — 1 = AW 2 , it will 
suffice to let \P/P - 1| < \W 2 \. 

Example: Let P = k/(s — 1), where k is in the range 2-5. We need to create a nominal 
model P = ko/(s — 1), with the smallest possible value of W 2 , which will not vary with 
frequency in this case. Two equations can be written using the above estimate, for the two 
extreme values of k, yielding k = 7/2, and W 2 = 3/7. 

For constructing the Nyquist plot, we observe that 

P(s)C(s) = (1 + A(s)W 2 (s))P(s)C(s). The path of the perturbed plant could be anywhere 
on a disk of radius \W 2 (s)P(s)C(s)\, centered on the nominal loci P(s)C(s). The robustness 
condition is that this disk should not intersect the critical point. This can be written as 



\1 + PC\ > \W 2 PC\ <— > 
\W 2 PC\ 

1 > TZ < 



1 + PC\ 

1 > \W 2 T\, (204) 

where T is the complementary sensitivity function. The last inequality is thus a condition 
for robust stability in the presence of multiplicative uncertainty parametrized with W 2 . 



18.5 Robust Performance 

The condition for good performance with plant uncertainty is a combination of the above two 
conditions. Graphically, the disk at the critical point, with radius \Wi\, should not intersect 
the disk of radius |W2-PC|, centered on the nominal locus PC. This is met if 

\WiS\ + \W 2 T\ < 1. (205) 

The robust performance requirement is related to the magnitude \PC\ at different frequen- 
cies, as follows: 

1. At low frequency, \WiS\ ~ \Wi/PC\, since \PC\ is large. This leads directly to the 
performance condition \PC\ > \Wi\ in this range. 

2. At high frequency, W 2 T\ ~ I^PC), since \PC\ is small. We must therefore have 
| PC | < l/\W 2 1, for robustness. 



18.6 Implications of Bode's Integral 
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18.6 Implications of Bode's Integral 

The loop transfer function PC cannot roll off too rapidly in the crossover region. The 
simple reason is that a steep slope induces a large phase loss, which in turn degrades the 
phase margin. To see this requires a short foray into Bode's integral. For a transfer function 
H(s), the crucial relation is 

1 Z" 00 d 

angle(H(ju )) = - / —(ln\H(ju)) ■ ln(coth(\u\/2))dv, (206) 

7T J-oo dV 

where v = ln(u/uo). The integral is hence taken over the log of a frequency normalized with 
Uq. It is not hard to see how the integral controls the angle: the function ln(coth(\v\/2)) is 
nonzero only near v — 0, implying that the angle depends only on the local slope d(ln\H\)/du. 
Thus, if the slope is large, the angle is large. 

Example: Suppose H(s) = uJq/s 71 , i.e., it is a simple function with n poles at the origin, 
and no zeros; cu is a fixed constant. It follows that \H\ = uJq/u 71 , and ln\H\ = —nln(u / ujq) , 
so that d{ln\H\)/dv = —n. Then we have just 

angle(H) = / ln(coth(\v\/2))dv = - — . (207) 

7T J-oo 2 

This integral is trivial to look up or compute. Each pole at the origin clearly induces 90° 
of phase loss. In the general case, each pole not at the origin induces 90° of phase loss for 
frequencies above the pole. Each zero at the origin adds 90° phase lead, while zeros not at 
the origin at 90° of phase lead for frequencies above the zero. In the immediate neighborhood 
of these poles and zeros, the phase may vary significantly with frequency. 
The Nyquist loci are clearly susceptible to these variations is phase, and the phase margin 
can be easily lost if the slope of PC at crossover (where the magnitude is unity) is too steep. 
The slope can safely be first-order (— 20dB /decade, equivalent to a single pole), and may be 
second-order 

(—4QdB I decade) if an adequate phase angle can be maintained near crossover. 

18.7 The Recipe for Loopshaping 

In the above analysis, we have extensively described what the open loop transfer function 
PC should look like, to meet robustness and performance specifications. We have said very 
little about how to get the compensator C, the critical component. For clarity, let the 
designed loop transfer function be renamed, L = PC. We will use concepts from optimal 
linear control for the MIMO case, but in the scalar case, it suffices to just pick 

C = L/P. (208) 
This extraordinarily simple step involves a plant inversion. 

The overall idea is to first shape L as a stable transfer function meeting the requirements of 
stability and robustness, and then divide through by the plant. 



92 



1 9 LINEAR Q UADRATIC REG ULATOR 



• When the plant is stable and has stable zeros (minimum-phase), the division can be 
made directly. 

• One caveat for the well-behaved plant is that lightly-damped poles or zeros should not 
be cancelled verbatim by the compensator, because the closed-loop response will be 
sensitive to any slight change in the resonant frequency. The usual procedure is to 
widen the notch or pole in the compensator, through a higher damping ratio. 

• Non-minimum phase or unstable behavior in the plant can usually be handled by 
performing the loopshaping for the closest stable model, and then explicitly considering 
the effects of adding the unstable parts. In the case of unstable zeros, we find that they 
impose an unavoidable frequency limit for the crossover. In general, the troublesome 
zeros must be faster than the closed-loop frequency response. 

In the case of unstable poles, the converse is true: The feedback system must be faster 
than the corresponding frequency of the unstable mode. 

When a control system involves multiple inputs and outputs, the ideas from scalar loop- 
shaping can be adapted using the singular value. We list below some basic properties of 
the singular value decomposition, which is analogous to an eigenvector, or modal, analysis. 
Useful properties and relations for the singular value are found in the section MATH FACTS. 
The condition for MIMO robust performance can be written in many ways, including a direct 
extension of our scalar condition 

W(W 1 S)+W(W 2 T) < 1. (209) 

The open-loop transfer matrix L should be shaped accordingly. In the following sections, we 
use the properties of optimal state estimation and control to perform the plant inversion for 
MIMO systems. 

19 LINEAR QUADRATIC REGULATOR 

19.1 Introduction 

The simple form of loopshaping in scalar systems does not extend directly to multivariable 
(MIMO) plants, which are characterized by transfer matrices instead of transfer functions. 
The notion of optimality is closely tied to MIMO control system design. Optimal controllers, 
i.e., controllers that are the best possible, according to some figure of merit, turn out to 
generate only stabilizing controllers for MIMO plants. In this sense, optimal control solutions 
provide an automated design procedure - we have only to decide what figure of merit to 
use. The linear quadratic regulator (LQR) is a well-known design technique that provides 
practical feedback gains. 



19.2 Full-State Feedback 
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19.2 Full-State Feedback 

For the derivation of the linear quadratic regulator, we assume the plant to be written in 
state-space form x = Ax + Bu, and that all of the n states x are available for the controller. 
The feedback gain is a matrix K, implemented as u = —K{x—x desired)- The system dynamics 
are then written as: 

X = (A - BK)X + BKXdesired- (210) 

Xdesired represents the vector of desired states, and serves as the external input to the closed- 
loop system. The "A- matrix" of the closed loop system is (A — BK), and the "B- matrix" 
of the closed-loop system is BK. The closed-loop system has exactly as many outputs as 
inputs: n. The column dimension of B equals the number of channels available in u, and 
must match the row dimension of K. Pole-placement is the process of placing the poles of 
(A — BK) in stable, suitably-damped locations in the complex plane. 

19.3 The Maximum Principle 

First we develop a general procedure for solving optimal control problems, using the calculus 
of variations. We begin with the statement of the problem for a fixed end time tf. 

choose u(t) to minimize J = tp(x(tf)) + / L(x(t),u(t),t)dt (211) 

J t 

subject to x = f(x(t),u(t),t) (212) 
x(t ) = x (213) 

where ip(x(tf),tf) is the terminal cost; the total cost J is a sum of the terminal cost and an 
integral along the way. We assume that L(x(t),u(t),t) is nonnegative. The first step is to 
augment the cost using the costate vector X(t). 



J = ip(x(t f ))+ [ tf (L + \ T (f-x))dt 

J t 



(214) 



Clearly, \(t) can be anything we choose, since it multiplies / — x — 0. Along the optimum 
trajectory, variations in J and hence J should vanish. This follows from the fact that J is 
chosen to be continuous in x, u, and t. We write the variation as 

SJ = ip x Sx(t f ) + / LJx + L u Su + X T f x Sx + X T f u Su - X T 5x dt, (215) 

where subscripts denote partial derivatives. The last term above can be evaluated using 
integration by parts as 



- f' X T 5xdt = -X T {t f )Sx(t f ) + X T {t )Sx(t ) + f tf X T Sxdt, (216) 

J t J t 



and thus 
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5J = ip x (x(t f ))5x(t f )+ [ tf (L u + \ T f u )5udt+ (217) 

J t 

f tf {L x + X T f x + X T )5xdt - X T {t f )Sx{t f ) + X T (t )8x{t ). 

J t 



Now the last term is zero, since we cannot vary the initial condition of the state by changing 
something later in time - it is a fixed value. This way of writing J makes it clear that there 
are three components of the variation that must independently be zero (since we can vary 
any of x, u, or x(tf)): 



L u + \ T f u = (218) 
L x + X T f x + X T = (219) 
^ x {x{t f ))-X T {t f ) = 0. (220) 

The second and third requirements are met by explicitly setting 

A T = -L x -X T f x (221) 
X T (t f ) = M<t f )). (222) 

The evolution of A is given in reverse time, from a final state to the initial. Hence we see the 
primary difficulty of solving optimal control problems: the state propogates forward in time, 
while the costate propogates backward. The state and costate are coordinated through the 
above equations. 

19.4 Gradient Method Solution for the General Case 

Numerical solutions to the general problem are iterative, and the simplest approach is the 
gradient method. It is outlined as follows: 

1. For a given x Q , pick a control history u(t). 

2. Propogate x = f(x,u,t) forward in time to create a state trajectory. 

3. Evaluate ip x (x(tf)), and the propogate the costate backward in time from tf to t a , 
using Equation 221. 

4. At each time step, choose Su = —K(L U + X T f u ), where K is a positive scalar or a 
positive definite matrix in the case of multiple input channels. 

5. Let u = u + 5u. 

6. Go back to step 2 and repeat loop until solution has converged. 



19.5 LQR Solution 
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The first three steps are consistent in the sense that x is computed directly from x(t Q and u, 
and A is computed directly from x and x(tf). All of S J in Equation 217 except the integral 
with Su is therefore eliminated explicitly. The choice of Su in step 4 then achieves SJ < 0, 
unless Su = 0, in which case the problem is solved. 

The gradient method is quite popular in applications and for complex problems may be the 
only way to effectively develop a solution. The major difficulties are computational expense, 
and the requirement of having a reasonable control trajectory to begin. 



19.5 LQR Solution 

In the case of the Linear Quadratic Regulator (with zero terminal cost), we set -0 = 0, and 

L = -x T Qx + ]-u T Ru, (223) 

where the requirement that L > implies that both Q and R are positive definite. In the 
case of linear plant dynamics also, we have 



fx 
fu 



x T Q 
u T R 
A 



(224) 
(225) 
(226) 
(227) 



so that 



x = Ax + Bu (228) 

x(t ) = x (229) 

A = -Qx - A T X (230) 

X(t f ) = (231) 

Ru + B T X = 0. (232) 

Since the systems are clearly linear, we try a connection A = Px. Inserting this into the A 
equation, and then using the x equation, and a substitution for u, we obtain 

PAx + A T Px + Qx- PBR- l B T Px + P = 0. (233) 

This has to hold for all x, so in fact it is a matrix equation, the matrix Riccati equation. 
The steady-state solution is given satisfies 



PA + A T P + Q- PBR' l B T P = 0. 



(234) 
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19.6 Optimal Full-State Feedback 

This equation is the matrix algebraic Riccati equation (MARE), whose solution P is needed 
to compute the optimal feedback gain K. The MARE is easily solved by standard numerical 
tools in linear algebra. 

The equation Ru + B T X = gives the feedback law: 

u = -R- 1 B T Px. (235) 



19.7 Properties and Use of the LQR 

Static Gain. The LQR generates a static gain matrix K, which is not a dynamical system. 
Hence, the order of the closed-loop system is the same as that of the plant. 

Robustness. The LQR achieves infinite gain margin: k g = oo, implying that the loci of 
(PC) (scalar case) or (det(I+PC) — 1) (MIMO case) approach the origin along the imaginary 
axis. The LQR also guarantees phase margin 7 > 60 degrees. This is in good agreement 
with the practical guidelines for control system design. 

Output Variables. In many cases, it is not the states x which are to be minimized, but 
the output variables y. In this case, we set the weighting matrix Q = C T Q'C, since y = Cx, 
and the auxiliary matrix Q' weights the plant output. 

Behavior of Closed-Loop Poles: Expensive Control. When R » C T Q'C, the cost 
function is dominated by the control effort u, and so the controller minimizes the control 
action itself. In the case of a completely stable plant, the gain will indeed go to zero, so that 
the closed-loop poles approach the open-loop plant poles in a manner consistent with the 
scalar root locus. 

The optimal control must always stabilize the closed-loop system, however, so there should 
be some account made for unstable plant poles. The expensive control solution puts stable 
closed-loop poles at the mirror images of the unstable plant poles. 

Behavior of Closed-Loop Poles: Cheap Control. When R « C T Q'C, the cost 
function is dominated by the output errors y, and there is no penalty for using large u. 
There are two groups of closed-loop poles. First, poles are placed at stable plant zeros, and 
at the mirror images of the unstable plant zeros. This part is akin to the high-gain limiting 
case of the root locus. The remaining poles assume a Butterworth pattern, whose radius 
increases to infinity as R becomes smaller and smaller. 

The Butterworth pattern refers to an arc in the stable left-half plane, as shown in the figure. 
The angular separation of n closed-loop poles on the arc is constant, and equal to 180°/n. 
An angle 90°/n separates the most lightly-damped poles from the imaginary axis. 



19.8 Proof of the Gain and Phase Margins 
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19.8 Proof of the Gain and Phase Margins 

The above stated properties of the LQR control are not difficult to prove. We need an 
intermediate tool (which is very useful in its own right) to prove it: the Liapunov function. 
Consider a scalar, continuous, positive definite function of the state V(x). If this function is 
always decreasing, except at the origin, then the origin is asymptotically stable. There is an 
intuitive feel to the Liapunov function; it is like the energy of a physical system. If it can 
be shown that energy is leaving such a system, then it will "settle." It suffices to find one 
Liapunov function to show that a system is stable. 
Now, in the case of the LQR control, pick 
V(x) = \x T Px. 

It can be shown that P is positive definite; suppose that instead of the design gain B, we 
have actual gain B 1 ', giving (with constant P based on the design system) 



V = 2x T Px (236) 
= 2x T P{Ax - B'R- 1 B T Kx) (237) 



x 



T < Q + PBR- 1 B T P)x-2x T KB'R- 1 B T Px. (238) 



where we used the Riccati equation in a substitution. Since we require V < for x ^ 0, we 
must have 

Q + PBG-2P(B- B')G> 0, (239) 

where G = R^ 1 B T P. Let B' = BN, where N is a diagonal matrix of random, unknown 
gains on the control channels. The above equation becomes 

Q - PB(I - 2N)R~ 1 B T P > 0. (240) 

This is satisfied if, for all channels, iV^j > 1/2. Thus, we see that the LQR provides for 
one-half gain reduction, and infinite gain amplification, in all channels. 
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The phase margin is seen (roughly) by allowing N to be a diagonal matrix of transfer 
functions N^i = In this case, we require that the real part of is greater than one- 
half. This is true for all |0j| < 60°, and hence sixty degrees of phase margin is provided in 
all channels. 



20 KALMAN FILTER 

20.1 Introduction 

In the previous section, we derived the linear quadratic regulator as an optimal solution 
for the full-state feedback control problem. The inherent assumption was that each state 
was known perfectly. In real applications, the measurements are subject to disturbances, 
and may not allow reconstruction of all the states. This state estimation is the task of a 
model-based estimator having the form: 

x = Ax + Bu + H(y - Cx) (241) 

The vector x represents the state estimate, whose evolution is governed by the nominal 
A and B matrices of the plant, and a correction term with the estimator gain matrix H. 
H operates on the estimation error mapped to the plant output y, since Cx = y. Given 
statistical properties of real plant disturbances and sensor noise, the Kalman Filter designs 
an optimal H. 




20.2 Problem Statement 

We consider the state-space plant model given by: 

x = Ax + Bu + Wx (242) 

y = Cx + W 2 . 

There are n states, m inputs, and I outputs, so that A has dimension n x n, B is n x m, 
and C is I x n. The plant subject to two random input signals, W\ and W<i- W\ represents 



20.3 Step 1: An Equation for t 
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disturbances to the plant, since it drives x directly; W 2 denotes sensor noise, which corrupts 
the measurement y. 

An important assumption of the Kalman Filter is that W\ and W 2 are each vectors of 
unbiased, independent white noise, and that all the n + I channels are uncorrelated. Hence, 
if E(-) denotes the expected value, 



E(W 1 (t 1 )W 1 (t 2 f) = V 1 8(t 1 -t 2 ) (243) 
E{W 2 (t 1 )W 2 {t 2 ) T ) = V^ih-h) (244) 
E{W 1 {t)W 2 {t) T ) = nxl . (245) 

Here 5(t) represents the impulse (or delta) function. V\ is an n x n diagonal matrix of 

intensities, and V 2 is an I x I diagonal matrix of intensities. 

The estimation error may be defined as e = x — x. It can then be verified that 

e = [Ax + Bu + Wx] - [Ax + Bu + H(y - Cx)} (246) 
= {A- HC)e + {W l - HW 2 ). 

The eigenvalues of the matrix A — HC thus determine the stability properties of the es- 
timation error dynamics. The second term above, W\ + HW 2 is considered an external 
input. 

The Kalman filter design provides H that minimizes the scalar cost function 

J = E(e T We), (247) 

where W is an unspecified symmetric, positive definite weighting matrix. A related matrix, 
the symmetric error covariance, is defined as 

E = £(ee T ). (248) 
There are two main steps for deriving the optimal gain H. 

20.3 Step 1: An Equation for t 

The evolution of S follows from some algebra and the convolution form of e(t). We begin 
with 



t = E{ee T + ee T ) (249) 
= E[(A - HC)ee T + (W 1 - HW 2 )e T + ee T (A T - C T H T ) + 
e(Wf-W?H T )}. 

The last term above can be expanded, using the property that 
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e(t) = e {A - HC)t e(0) + f e^-^^HwAr) - HW 2 (r)) dr. 

Jo 

We have 



E{e{W? - W?H T )) = e( A - HC »E(e(0)(W? -W?H T )) + 

J 1 e ^-HC)(t-r) E ^ Wi ^ _ Mj ( T )j 

{W T {t) _ wZ(t)H T )) dr 

= f e (A-HC)(t-r) E U Wl ( T ) _ HW 2 (T)) 

Jo 

(WTq _ W ?(t)H T )) dr 

= C e (A-HC)(t-r) ,y iS , t _ r) + H y 2H T 5 u — t\) dr 
Jo 

1 1 rp 

= - Vl + -HV 2 H T . 

To get from the first right-hand side to the second, we note that the initial condition e(0) is 
uncorrected with W[ — W^H 7 . The fact that W 1 and HW 2 are uncorrelated leads to the 
third line, and the final result follows from 



i.e., the written integral includes only half of the impulse. 

The final expression for E(e(W± — W^H 7 )) is symmetric, and therefore appears in Equa- 
tion 249 twice, leading to 

t = (A- HC)E + Tj(A t - C T H T ) + Vx + HV 2 H T . (250) 

This equation governs propagation of the error covariance. It is independent of the initial 
condition e(0), and depends on the (as yet) unknown estimator gain matrix H. 

20.4 Step 2: H as a Function of E 

We now make the connection between S = E(ee T ) (a matrix) and J = E(e T We) (a scalar). 
The trace of a matrix is the sum of its diagonal elements, and it can be verified that 

J = E(e T We) = trace(EW). (251) 
We now introduce an auxiliary cost function defined as 

J' = trace(EW + AF) , (252) 



20.5 Properties of the Solution 
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where F is an n x n matrix of zeros, and A is an n x n matrix of unknown Lagrange multipliers. 
Note that since F is zero, J' = J, so minimizing J' solves the same problem. Lagrange 
multipliers provide an ingenious mechanism for drawing constraints into the optimization; 
the constraint we invoke is the evolution of S, Equation 250: 

J' = trace (ZW + A(-S + AS - HCH + SA T - Y>C T H T + V 1 + HV 2 H T )) (253) 

If J' is an optimal cost, it follows that dJ'/dH = 0, i.e., the correct choice of H achieves 
an extremal value. We need the following lemmas, which give the derivative of a trace with 
respect to a constituent matrix: 



-^-trace(-AHCE) = -A T SC T 
oH 



A trace( _ ACT) = _ AEC - 

3H 

d 

—trace(AHV 2 H T ) = A T HV 2 + AHV 2 . 
on 

Proofs of the first two are given at the end of this section; the last lemma uses the chain 
rule, and the previous two lemmas. Next, we enforce A = A T , since the values are arbitrary. 
Then the condition dJ' /dH = leads to 



= 2A(-SC T + HV 2 ), satisfied if 

H = EC T V 2 -\ (254) 

Hence the estimator gain matrix H can be written as a function of S. Inserting this back 
into Equation 250, we obtain 

t = AS + SA T + V 1 - EC r V r 2 _1 CE. (255) 

Equations 254 and 255 represent the practical solution to the Kalman filtering problem, 
which minimizes the squared-norm of the estimation error. The evolution of S is always 
stable, and depends only on the constant matrices [A, C, Vi, V 2 \. Notice also that the result 
is independent of the weighting matrix W, which might as well be the identity. 



20.5 Properties of the Solution 

The solution involves a matrix Riccati equation, like the LQR, suggesting a duality with the 
LQR problem. This is in fact the case, and the same analysis and numerical tools can be 
applied to both methodologies. 

The steady-state solution for S is valid for time-invariant systems, leading to a more common 
MARE form of Equation 255: 



= AS + SA T + V 1 - Y>C T V 2 l CY>. 



(256) 
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Duality of Linear Quadratic Regulator and Kalman Filter 



Linear Quadratic Regulator 


Kalman Filter 


x = Ax + Bu 
u = —Kx 


x = Ax + Bu + Wi 

y = Cx + W 2 

x = Ax + Bu + H(y - Cx) 


2J = J™(x T Qx + u r Ru)dt 
Q > 0,R > 


J = E(e r We) 
Vi > 0, V 2 > 


K = R- l B T P 

PA + A T P + Q- PBR~ X B T P = 


H = EC T V 2 - L 

TjA t + AS + V x - EC^V^CE = 



The Kalman Filter is guaranteed to create a stable nominal dynamics A — HC, as long as 
the plant is fully state-observable. This is dual to the stability guarantee of the LQR loop, 
when the plant is state-controllable. Furthermore, like the LQR, the KF loop achieves 60° 
phase margin, and infinite gain margin, for all the channels together or independently. 
The qualitative dependence of the estimator gain H = TiC T V 2 ~ 1 on the other parameters 
can be easily seen. Recall that V\ is the intensity matrix of the plant disturbance, V 2 is the 
intensity of the sensor noise, and S is the error covariance matrix. 

• A large uncertainty S creates large H, placing emphasis on the corrective action of the 
filter. 

• A small disturbance Vi, and large sensor noise V 2 creates a small H, weighting the 
model dynamics Ax + Bu more. 

• A large disturbance V±, and small sensor noise V 2 creates a large H, so that the filter's 
correction is dominant. 

The limiting closed-loop poles of the Kalman filter are similar, and dual to those of the LQR: 

• V 2 « Vi 1 . good sensors, large disturbance, H » 1, dual to cheap-control problem. 
Some closed-loop poles go to the stable plant zeros, or the mirror image of unstable 
plant zeros. The remaining poles follow a Butterworth pattern whose radius increases 
with increasing Vi/V 2 . 

• V 2 » Vi. poor sensors, small disturbance, if small, dual to expensive-control problem. 
Closed-loop poles go to the stable plant poles, and the mirror images of the unstable 
plant poles. 

20.6 Combination of LQR and KF 

An optimal output feedback controller is created through the use of a Kalman filter coupled 
with an LQR full-state feedback gain. This combination is usually known as the Linear 
Quadratic Gaussian design, or LQG. For the plant given as 



20. 7 Proofs of the Intermediate Results 
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x = Ax + Bu + Wi 

y = CX + W 2 , 

we put the Kalman Filter and controller gain G together as follows: 



k = Ax + Bu + H(y-Cx) (257) 
u = -Kx. (258) 



B 



C(s) 



-K 



()) = (sI-A)" 1 



There are two central points to this construction: 

1. Separation Principle: The eigenvalues of the nominal closed-loop system are made 
of up the eigenvalues of (A — HC) and the eigenvalues of (A — BK), separately. See 
proof below. 

2. Output Tracking: This compensator is a stand-alone system that, as written, tries to 
drive its input y to zero. It can be hooked up to receive tracking error e(s) = r(s) — y(s) 
as an input instead, so that it is not limited to the regulation problem alone. In this 
case, x no longer represents an estimated state, but rather an estimated state tracking 
error. We use the output error as a control input in the next section, on loopshaping 
via loop transfer recovery. 

20.7 Proofs of the Intermediate Results 
20.7.1 Proof that E(e T We) = trace(Y<W) 



E(e T We) = B Ee, D% 
\i=i \j=i 

n n 

= EEw 

i=l j=\ 
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the transpose of W being valid since it is symmetric. Now consider the diagonal elements of 
the product T,W: 



EnWii + £ 12 W 21 + • • • 



Z 21 W l2 + Z 22 W 22 + ■■■ ■ 



trace(EW) = EE S u^- D 
i=i 3=1 

20.7.2 Proof that ^trace(-AHCT,) = -A T EC q 



trace(AHB) = trace 



n I 

E A? H H jk B k i 

j=l k=l 
n n I 

EE^'J H H jk B ki , 
i=lj=l k=l 



, the i/'th element 



where the second form is a sum over % of the zz'th elements. Now 



d 



OH, 



d 



-trace(AHB) 



—trace(AHB) 
oH 



— A ijo B koi 

i=l 

= (BA) kojo 
= {BA) T joko ^ 

= {BAY 

= A T B T . □ 



20.7.3 Proof that ^trace{-AEC T H T ) = -AEC 1 



trace(AH T ) = trace 



n 



= trace 



3=1 

n 



, the il'th. element 



E A i3 H U 



— E E ^ijBij j 
i=l j=l 



where the last form is a sum over the iz'th elements. It follows that 
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trace(AH T ) = A io]o 



dH iojo 



d 

—trace(AH T ) = A. □ 
oH 

20.7.4 Proof of the Separation Principle 

Without external inputs W\ and W2, the closed-loop system evolves according to 



d 



f x ] 




I * J 





A -BK 
HC A-BK-HC 




\ x ) 




\ - J 






Using the definition of estimation error e = x — x, we can write the above in another form: 



A- BK BK 
A- EC 

If A' represents this compound A- matrix, then its eigenvalues are the roots of det(sl — A') = 
0. However, the determinant of an upper triangular block matrix is the product of the 
determinants of each block on the diagonal: det(sl — A') = det(sl — (A — BK))det(sI — 
(A — HC)), and hence the separation of eigenvalues follows. 

21 LOOP TRANSFER RECOVERY 

21.1 Introduction 

The Linear Quadratic Regulator(LQR) and Kalman Filter (KF) provide practical solutions 
to the full-state feedback and state estimation problems, respectively. If the sensor noise and 
disturbance properties of the plant are indeed well-known, then an LQG design approach, 
that is, combining the LQR and KF into an output feedback compensator, may yield good 
results. The LQR tuning matrices Q and R would be picked heuristically to give a reasonable 
closed-loop response. 

There are two reasons to avoid this kind of direct LQG design procedure, however. First, 
although the LQR and KF each possess good robustness properties, there do exist plants 
for which there is no robustness guarantee for an LQG compensator. Even if one could 
steer clear of such pathological second problem is that this design technique has no 

clear equivalent in frequency space. It cannot be directly mapped to the intuitive ideas of 
loopshaping and the Nyquist plot, which are at the root of feedback control. 
We now reconsider just the feedback loop of the Kalman filter. The KF has open-loop 
transfer function L(s) = C(p(s)H, where (p(s) = (si — A)^ 1 . This follows from the estimator 
evolution equation 
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O 



H 



(|)(S) 



C 



'x = Ax + Bu + H(y - Cx) 

and the figure. Note that we have not included the factor Bu as part of the figure, since it 
does not affect the error dynamics of the filter. 

As noted previously, the KF loop has good robustness properties, specifically to perturbations 
at the output y, and further is amenable to output tracking. In short, the KF loop is an 
ideal candidate for a loopshaping design. Supposing that we have an estimator gain H which 
creates an attractive loop function L(s), we would like to find the compensator C(s) that 
establishes 



P(s)C(s) « C<p(s)H,oi (259) 
C<j){s)BC{s) « C<j>(s)H. 

It will turn out that the LQR can be set up so that the an LQG-type compensator achieves 
exactly this result. The procedure is termed Loop Transfer Recovery (LTR), and has two 
main parts. First, one carries out a KF design for H, so that the Kalman filter loop itself has 
good performance and robustness properties. In this regard, the KF loop has sensitivity func- 
tion S(s) = (I + C(j)(s)H)^ 1 and complementary sensitivity T(s) = (I + C(j)(s)H)^ 1 C(f)(s)H . 
The condition a(Wi(s) S (s)) + W(W 2 (s)T(s)) < 1 is sufficient for robust performance with 
multiplicative plant uncertainty at the output. Secondly, we pick suitable parameters of the 
LQR design, so that the LQG compensator satisfies the approximation of Equation 259. 
LTR is useful as a SISO control technique, but has a much larger role in multivariable control. 

21.2 A Special Property of the LQR Solution 

Letting Q = C T C and R = pi, where / is the identity matrix, we will show (roughly) that 



Hm(Vptf) = WC, 

where K is the LQR gain matrix, and W is an orthonormal matrix, for which W T W = I. 
First recall the gain and Riccati equations for the LQR: 



K = R- l B T P 

= Q + PA + A T P - PBR- l B T P. 



21.3 The Loop Transfer Recovery Result 
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Now Q = C T C = C T W T WC = (WC) T WC. The Riccati equation becomes 

= p(WC) T WC + pPA + P A T P - PBB T P = 0. 
In the limit as p — > 0, it must be the case that P — > also, and so in this limit 

p(WC) T WC « PBB T P 

= (B T P) T B T P 

= (R- 1 B T P) T RR(R- 1 B T P) 

= p 2 K T K — > 

WC « y/pK. □ 

Note that another orthonormal matrix W could be used in separating K T from K in the 
last line. This matrix may be absorbed into W through a matrix inverse, however, and so 
does not need to be written. The result of the last line establishes that the plant must be 
square: the number of inputs (i.e., rows of K) is equal to the number of outputs (i.e., rows 
ofC). 

Finally, we note that the above property is true only for LQR designs with minimum-phase 
plants, i.e., those with only stable zeros (Kwakernaak and Sivan). 

21.3 The Loop Transfer Recovery Result 

The theorem is stated as: If \im p ^ (y/pK) = WC (the above result), with W an orthonormal 
matrix, then the limiting LQG controller C(s) satisfies 

lim P(s)C(s) = C<P{s)H. 

The LTR method is limited by two conditions: 

• The plant has an equal number of inputs and outputs. 

• The design plant has no unstable zeros. The LTR method can be in fact be applied in 
the presence of unstable plant zeros, but the recovery is not to the Kalman filter loop 
transfer function. Instead, the recovered function will exhibit reasonable limitations 
inherent to unstable zeros. See Athans for more details and references on this topic. 

The proof of the LTR result depends on some easy lemmas, given at the end of this section. 
First, we develop C(s), with the definitions <f>(s) = (si- A)- 1 andX(s) = ((j>- 1 (s)+HC)- 1 = 
(si - A + HC)- 1 . 
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C(s) = K(sl - A + BK + HC)- 1 H 

= K(X^ 1 (s) + BK)~ 1 H, then use Lemma 2 — >■ 

= K(X(s) - X(s)B(I + KX{s)BY 1 KX(s))H 

= KX{s)H - KX(s)B(I + KX(s)BY l KX(s)H 

= (I- KX(s)B(I + KX(s)B)- 1 )KX(s)H, then use Lemma 3 ->■ 

= (/ + KX(s)B)- 1 KX(s)H 

= (^pI+^pKX(s)B)' 1 ^pKX(s)H. 

Next we invoke the result from the LQR design, with p — > 0, to eliminate y/pK: 



lim C{s) = (WCX(s)B)- 1 WCX(s)H 
= {CX{s)B)- 1 CX{s)H. 

In the last expression, we used the assumption that W is square and invertible, both prop- 
erties of orthonormal matrices. Now we look at the product CX(s): 



CX(s) = C{SI - A + HC)- 1 

= C(0 _1 (s) + HO)' 1 , then use Lemma 2 -»■ 

= C{(j){s)- (j)(s)H(I + C(j){s)H)- 1 C(j){s)) 

= (I- C<j)(s)H{I + C0(s)i/)~ 1 )C0(s), then use Lemma 3 -> 

= (I + C^H)- 1 ^). 

This result, reintroduced into the limiting compensator, gives 



limC(s) = ((/ + C0(s) J F/)~ 1 C^(s) J B)- 1 (/ + C0(s)F)^ 1 C^(s)i/ 

= (C0(s) J B)- 1 C(^(s)i/ 
= p- 1 {s)C(j>{s)H. 

Finally it follows that \im p ^ P(s)C(s) = C(f>(s)H, as desired. 



21.4 Usage of the Loop Transfer Recovery 

The idea of LTR is to "recover" a Kalman filter loop transfer function L(s) = C<p(s)H, by 
using the limiting cheap-control LQR design, with Q = C T C and R = pi. The LQR design 
step is thus trivial. 
Some specific techniques are useful. 



21.5 



Three Lemmas 
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• Scale the plant outputs (and references), so that one unit of error in one channel is as 
undesirable as one unit of error in another channel. For example, in depth and pitch 
control of a large submarine, one meter of depth error cannot be compared directly 
with one radian of pitch error. 

• Scale the plant inputs in the same way. One Newton of propeller thrust cannot be 
compared with one radian of rudder angle. 

• Design for crossover frequency. The bandwidth of the controller is roughly equal to 
the frequency at which the (recovered) loop transfer function crosses over OdB. Often, 
the bandwidth of is a more intuitive design parameter than is, for example, the high- 
frequency multiplicative weighting W2- Quantitative uncertainty models are usually at 
the cost of a lengthy identification effort. 

• Integrators should be part of the KF loop transfer function, if no steady-state error 
is to be allowed. Since the Kalman filter loop has only as many poles as the plant, 
the plant input channels must be augmented with the necessary additional poles (at 
the origin). Then, once the KF design is completed, and the compensator C(s) is 
constructed, the integrators are moved from the plant over to the input side of the 
compensator. The tracking errors will accrue as desired. 

21.5 Three Lemmas 
Lemma 1: Matrix Inversion 



(A + BCD) 



-i 



A 



-i 



A^BiC- 1 + DA- 1 B)- 1 DA- 1 . 



Proof: 



(A + BCD)(A + BCD) 
A(A + BCD) 
{A + BCD) 



i 



I 

I - BCD{A + BCD)' 1 

A' 1 - A- l BCD(A + BCD)' 1 

A' 1 - A~ X BCD{1 + A' 1 BCD)' 1 A' 1 

A' 1 - A^B^D^C- 1 + A- 1 B)- 1 A- 1 

A' 1 - A^B^C- 1 + DA- 1 B)- 1 DA- 1 . □ 



i 



i 



Lemma 2: Short Form of Lemma 1 



(X- 1 + BD) 



-i 



X - XB(I + DXB)- 1 DX 



Proof: substitute A = X 1 and C = I into Lemma 1. 



Lemma 3 
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I-A(I + A) 



-i 



-i 



Proof: 



I-A(I + A) 



-i 



(/ + A){I + A)- 1 - 
(I + A-A){I + A) 
(I + A)- 1 . □ 



A(I + A) 



-i 



-i 
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22.1 Vectors 



22.1.1 Definition 



A vector has a dual definition: It is a segment of a a line with direction, or it consists of its 
projection on a reference system Oxyz, usually orthogonal and right handed. The first form 
is independent of any reference system, whereas the second (in terms of its components) 
depends directly on the coordinate system. Here we use the second notation, i.e., x is meant 
as a column vector, whose components are found as projections of an (invariant) directed 
segment on a specific reference system. 

We use the overhead arrow to denote a column vector, i.e., a linear segment with a direction. 
For example, in three-space, we write a vector in terms of its components with respect to a 
reference system as 



The elements of a vector have a graphical interpretation, which is particularly easy to see in 
two or three dimensions. 




z 



7 



a 



22.1 Vectors 
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1. Vector addition: 



r 2 1 

i 

7 



> + 



a + b 

3 1 

3 
2 



r 5 1 

4 
9 



Graphically, addition is stringing the vectors together head to tail. 
2. Scalar multiplication: 









2 x < 


[il- 








{ "14 J 



22.1.2 Vector Magnitude 

The total length of a vector of dimension m, its Euclidean norm, is given by 



\x\ 



This scalar is commonly used to normalize a vector to length one. 
22.1.3 Vector Dot or Inner Product 

The dot product of two vectors is a scalar equal to the sum of the products of the corre- 
sponding components: 



x-y = x T y = Y, x iVi- 

i=l 



The dot product also satisfies 



x-y=\\x\\\\y\\cose, 
where 6 is the angle between the vectors. 
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22.1.4 Vector Cross Product 

The cross product of two three-dimensional vectors x and y is another vector z, x x y = z, 
whose 

1. direction is normal to the plane formed by the other two vectors, 

2. direction is given by the right-hand rule, rotating from x to y, 

3. magnitude is the area of the parallelogram formed by the two vectors - the cross 
product of two parallel vectors is zero - and 

4. (signed) magnitude is equal to sin#, where 9 is the angle between the two 
vectors, measured from x to y. 

In terms of their components, 



22.2 Matrices 
22.2.1 Definition 

A matrix, or array, is equivalent to a set of column vectors of the same dimension, arranged 
side by side, say 



This matrix has three rows (m = 3) and two columns (n — 2); a vector is a special case of a 
matrix with one column. Matrices, like vectors, permit addition and scalar multiplication. 
We usually use an upper-case symbol to denote a matrix. 

22.2.2 Multiplying a Vector by a Matrix 

If denotes the element of matrix A in the z'th row and the j'th column, then the multi- 
plication c = Av is constructed as: 



x x y = x\ x 2 
Vi V2 




(x 2 yz - x 3 y 2 )i 
(x 3 y! - X!y 3 )j > 
(x x y 2 - x 2 yi)k 



A — [a b] 



2 3 
1 3 
7 2 



n 



d = A a vi + A i2 v 2 H h A in v, 



12 A a v i 



where n is the number of columns in A. c will have as many rows as A has rows (m). Note 
that this multiplication is defined only if v has as many rows as A has columns; they have 
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consistent inner dimension n. The product vA would be well-posed only if A had one row, 
and the proper number of columns. There is another important interpretation of this vector 
multiplication: Let the subscript : indicate all rows, so that each A : j is the j'th column 
vector. Then 

c — Av = A. l v l + A. 2 v 2 H h A. n v n . 

We are multiplying column vectors of A by the scalar elements of v. 

22.2.3 Multiplying a Matrix by a Matrix 

The multiplication C = AB is equivalent to a side-by-side arrangement of column vectors 
C : j = AB.j, so that 

C = AB = [AB :1 AB. 2 ■■■ AB. k ], 

where k is the number of columns in matrix B. The same inner dimension condition applies 
as noted above: the number of columns in A must equal the number of rows in B. Matrix 
multiplication is: 

1. Associative. (AB)C = A(BC). 

2. Distributive. A(B + C) = AB + AC, (B + C)A = BA + CA. 

3. NOT Commutative. AB ^ BA, except in special cases. 



22.2.4 Common Matrices 

Identity. The identity matrix is usually denoted /, and comprises a square matrix with 
ones on the diagonal, and zeros elsewhere, e.g., 



1 











1 











1 



The identity always satisfies AI nxn = I mxm A = A. 



Diagonal Matrices. A diagonal matrix is square, and has all zeros off the diagonal. For 
instance, the following is a diagonal matrix: 



4 











-2 











3 



The product of a diagonal matrix with another diagonal matrix is diagonal, and in this case 
the operation is commutative. 
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22.2.5 Transpose 

The transpose of a vector or matrix, indicated by a T superscript results from simply swap- 
ping the row-column indices of each entry; it is equivalent to "flipping" the vector or matrix 
around the diagonal line. For example, 



A 




{1 2 3} 



14 8 
2 5 9 



A very useful property of the transpose is 



(AB) T = B T A T . 

22.2.6 Determinant 

The determinant of a square matrix A is a scalar equal to the volume of the parallelepiped 
enclosed by the constituent vectors. The two-dimensional case is particularly easy to re- 
member, and illustrates the principle of volume: 



det 



det(A) = A U A 2 2 - A 2 iA 12 
1 -1 



1 1 



1 + 1 = 2. 









\ Area = det(A) 




'/:■ . 



In higher dimensions, the determinant is more complicated to compute. The general formula 
allows one to pick a row k, perhaps the one containing the most zeros, and apply 



j=n 

det(A) = Y,A kj (-l) k+ >A kj , 
i=i 

where A k j is the determinant of the sub-matrix formed by neglecting the fc'th row and the 
j'th column. The formula is symmetric, in the sense that one could also target the /c'th 
column: 
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det(A) = Y,A jk (-l) k+ >A jk . 
i=i 

If the determinant of a matrix is zero, then the matrix is said to be singular - there is no 
volume, and this results from the fact that the constituent vectors do not span the matrix 
dimension. For instance, in two dimensions, a singular matrix has the vectors colinear; in 
three dimensions, a singular matrix has all its vectors lying in a (two-dimensional) plane. 
Note also that det(A) = det(A T ). If det(A) ^ 0, then the matrix is said to be nonsingular. 

22.2.7 Inverse 

The inverse of a square matrix A, denoted A~ 1 , satisfies AA' 1 = A~ 1 A = I. Its computation 
requires the determinant above, and the following definition of the n x n adjoint matrix: 



adj(A) 



(-l) 1+1 A n ••• (-l) 1+ "A ln 

(-l)" +1 A nl ••• (-l)" + "A nn . 
Once this computation is made, the inverse follows from 



-i T 



A' 



adj(A) 
det(A) 



If A is singular, i.e., det(A) = 0, then the inverse does not exist. The inverse finds common 
application in solving systems of linear equations such as 



Ax = b — > x = A~ l b. 



22.2.8 Trace 

The trace of a square matrix is the sum of the diagonals: 

n 

tr(A) = Y,An. 
i=i 

22.2.9 Eigenvalues and Eigenvectors 

A typical eigenvalue problem is stated as 



Ax = Ax, 
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where A is an n x n matrix, x is a column vector with n elements, and A is a scalar. We ask 
for what nonzero vectors x (right eigenvectors), and scalars A (eigenvalues) will the equation 
be satisfied. Since the above is equivalent to (A — XI)x = 0, it is clear that det(A — XI) = 0. 
This observation leads to the solutions for A; here is an example for the two-dimensional 

case: 



A - XI 
det(A - XI) 



4 
2 



-5 
-3 



4 - A -5 
2 -3- A 

(4- A)(-3- A) + 10 
A 2 - A - 2 
(A + l)(A-2). 



Thus, A has two eigenvalues, Ai 
x. In this example, 



-1 and A2 = 2. Each is associated with a right eigenvector 



(A - X 1 I)x 1 

' 5 -5 
2 -2 



Xi 





— * 





Xi = {V2/2, V2/2} J 



(A 
' 2 
2 



- X 2 I)x 2 
-5 



-5 



£2 





— > 





x 2 = {5^/29, 2^/29} 



Eigenvectors are defined only within an arbitrary constant, i.e., if x is an eigenvector then cx 
is also an eigenvector for any c 7^ 0. They are often normalized to have unity magnitude, and 
positive first element (as above). The condition that rank(A — Xil) = rank(A) — 1 indicates 
that there is only one eigenvector for the eigenvalue A^; more precisely, a unique direction 
for the eigenvector, since the magnitude can be arbitrary. If the left-hand side rank is less 
than this, then there are multiple eigenvectors that go with Aj. 

The above discussion relates only the right eigenvectors, generated from the equation Ax = 
Ax. Left eigenvectors, defined as if 1 A = Xy 1 ^, are also useful for many problems, and can 
be defined simply as the right eigenvectors of A T . A and A T share the same eigenvalues A, 
since they share the same determinant. Example: 



(A T -X 1 I)y 1 = 
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(A T -X 2 I)y 2 = 




T 



22.2.10 Modal Decomposition 

For simplicity, we consider matrices that have unique eigenvectors for each eigenvalue. The 
right and left eigenvectors corresponding to a particular eigenvalue A can be defined to have 
unity dot product, that is xjyi = 1, with the normalization noted above. The dot products 
of a left eigenvector with the right eigenvectors corresponding to different eigenvalues are 
zero. Thus, if the set of right and left eigenvectors, V and W, respectively, is 

V = [xi • • • x n ] , and 
W = [yi---y n ], 

then we have 

W T V = I, or 

w T = v- 1 . 

Next, construct a diagonal matrix containing the eigenvalues: 



A = 



Ai 








A n 



it follows that 

AV = VA — > 
A = VAW T 

n 

= ^A^wf. 

i=l 

Hence A can be written as a sum of modal components. 2 

2 By carrying out successive multiplications, it can be shown that A k has its eigenvalues at A^, and keeps 
the same eigenvectors as A. 
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22.2.11 Singular Value 

Let G be an mx n real or complex matrix. The singular value decomposition (SVD) computes 
three matrices satisfying 



G = UEV% 

where U is m x m, Eismxn, and V is nxn. The star notation indicates a complex-conjugate 
transpose (the Hermitian of the matrix) . The matrix E has the form 



Si 




where Ei 



(71 

• 

a, 



v J 



and p 



minim, n) 



Each nonzero entry on the diagonal of matrix E x is a real, positive 



singular value, ordered such that (j\ > a 2 > • • • a p . Each of is an eigenvalue of G H G (or 
GG H ). The notation is common that <j\ = a, the maximum singular value, and a p = a, 
the minimum singular value. The auxiliary matrices U and V are unitary, i.e., they satisfy 
X* = X^ 1 . They are defined explicitly: U is the matrix of right eigenvectors of GG H , and 
V is the matrix of right eigenvectors of G H G. Like eigenvalues, the singular values of G 
are related to projections. Oi represents the Euclidean size of the matrix G along the z'th 
singular vector: 



a = max\\ x \\=i\\Gx 
g_ = min\\ x \\ = i\ \Gx\ 

Other properties of the singular value include: 

• W(AB) < a(A)a(B). 

• a(A) = ^X max (A*A). 

• a(A) = ^X min (A*A). 

• a(A) = IMA- 1 ). 

• a{A) = 1/aiA- 1 ). 

22.3 Laplace Transform 
22.3.1 Definition 



The Laplace transform projects time-domain signals into a complex frequency- domain equiv- 
alent. The signal y(t) has transform Y(s) defined as follows: 



22.3 Laplace Transform 
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poo 

Y(s) = L(y(t)) = / y(r)e-^dr, 
Jo 



where s is a complex variable, properly constrained within a region so that the integral 
converges. Y(s) is a complex function as a result. Note that the Laplace transform is linear, 
and so it is distributive: L(x(t) + y(t)) = L(x(t)) + L(y(t)). The following table gives a list 
of some useful transform pairs and other properties, for reference. 

The last two properties are of special importance: for control system design, the differenti- 
ation of a signal is equivalent to multiplication of its Laplace transform by s; integration of 
a signal is equivalent to division by s. The other terms that arise will cancel if y(0) = 0, or 
if y(0) is finite. 

22.3.2 Convergence 

We note first that the value of s affects the convergence of the integral. For instance, if 
y(t) = e*, then the integral converges only for Re(s) > 1, since the integrand is e l ~ s in this 
case. Although the integral converges within a well-defined region in the complex plane, the 
function Y(s) is defined for all s through analytic continuation. This result from complex 
analysis holds that if two complex functions are equal on some arc (or line) in the complex 
plane, then they are equivalent everywhere. It should be noted however, that the Laplace 
transform is defined only within the region of convergence. 

22.3.3 Convolution Theorem 

One of the main points of the Laplace transform is the ease of dealing with dynamic systems. 
As with the Fourier transform, the convolution of two signals in the time domain corresponds 
with the multiplication of signals in the frequency domain. Consider a system whose impulse 
response is g(t), being driven by an input signal x(t); the output is y(t) = g(t) * x(t). The 
Convolution Theorem is 




Here's the proof given by Siebert: 



Y(s) = f 
Jo 



CO 



y{t)e~ st dt 




x(t) / g(t - t) h(t - t) e- st dt dr 



/ g(t - t) h(t - t) x(t) dr e~ st dt 
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y(t) <r 

(Impulse) 5(t) <- 

(Unit Step) 1(f) <- 

(Unit Ramp) t -f- 



-at 



b — a 



1 

ab 



1 + 



1 



(be- at - ae~ bt ) 



OJ r , 



a — b 



sin uj, 



1 - 



sin [uj n Jl- ( 2 t + (p) i 



4> = tan" 



c 



sin u;i -f- 



cos cut -f- 



e Q * sin ut <- 



cos u;t <r 



-e- ht ) <r 



(Time Derivative) 



dt 



(Time Integral) / y(r)d 
Jo 



t i — y 



Y(s) 

1 
1 

s 
1 



s + a 

UJ 

S 2 + UJ 2 

s 

S 2 + UJ 2 
UJ 



(s + a) 2 + uj 2 
s + a 

(s + a) 2 + uj 2 
1 

(s + a)(s + b) 
1 

s(s + a)(s + b) 

s 2 + 2(u n s + 

s(s 2 + 2Co; n s + uj 2 ) 



(Pure Delay) y(t - r)l(t - r) « — ► Y(s)e" 

efy(t) 



^— ► s y( 5 )-y(o) 



y(g) | gly(f}dt 



22.4 Background for the Mapping Theorem 
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J o 



roo 

/ x(t) G(s)e- ST dr 



G(s)X(s), 



where h{t) is the unit step function. When g{t) is the impulse response of a dynamic system, 
then y{t) represents the output of this system when it is driven by the external signal x(t). 

22.3.4 Solution of Differential Equations by Laplace Transform 

The Convolution Theorem allows one to solve (linear time-invariant) differential equations 
in the following way: 

1. Transform the system impulse response g(t) into G(s), and the input signal x(t) into 
X(s), using the transform pairs. 

2. Perform the multiplication in the Laplace domain to find Y(s). 

3. Ignoring the effects of pure time delays, break Y (s) into partial fractions with no powers 
of s greater than 2 in the denominator. 

4. Generate the time-domain response from the simple transform pairs. Apply time delay 
as necessary. 

Specific examples of this procedure are given in a later section on transfer functions. 
22.4 Background for the Mapping Theorem 

The mapping theorem uses concepts from complex analysis, specifically Cauchy's Theorem 
and the Residue Theorem. References for this section include Ogata and Hildebrand. 
First, consider a function of the complex variable s = u + iv: f(s). We say f(s) is analytic in 
a region S if it has finite derivative and takes only one value at each point s in S. Therefore 
discontinuous or multi- valued functions, e.g., y/s, are not analytic functions. Polynomials in 
s are analytic, as are many functions that can be written as a Taylor or Laurent series. An 
especially important class for control system design is the rational function, a polynomial in 
s divided by another polynomial in s. Rational functions are consequently zero for certain 
values of s, the roots of the numerator, and undefined for other values of s, the roots of the 
numerator, also called the poles. 
The integral of interest here is 



taken on a path in the s-plane. A closed path in the complex s-plane leads to a closed path 
in the f(s) plane, but more than one point in the s plane can map to a single /(s)-plane 
point, so the number of complete loops may not be the same. 

The usual rules of integration apply in complex analysis, so that, insofar as the antiderivative 
of f(s), denoted F(s) exists, and f(s) is analytic on the path, we have 
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f S2 f(s)ds = F(s 2 )-F( Sl ). 

J Sl 

It appears that this integral is zero for a closed path, since si = S2- Indeed, Cauchy's theorem 
states that it is, provided that f(s) is analytic on the path, and everywhere within the region 
enclosed. This latter condition results from the following observation. Consider the function 
f(s) = as n ; on a circular path of radius r, we have s = re td , and thus 

iar n+1 r e^ n+1 ^d9 
Jo 

iar n+1 [ n \cos(n + 1)0 + % sin(n + l)9}d9. 
Jo 

The second integral is clearly zero for all n, whereas the first is zero except in the case of 
n = — 1, for which we obtain 



Jf(s)ds = 



J f(s)ds = ia2n. 

This result does not depend on r, and so applies to a vanishingly small circle around the 
point s = 0. It can be shown also that the result holds for any closed path around the simple 
pole at s = 0, which characterizes the function. The residue theorem is an extension to an 
arbitrary number of simple poles which are enclosed by a path: 

J f(s)ds = i2irJ2 a i' 
The constants a { are the residues associated with the poles, i.e., 



/W = _55_ + _5S_ + ... 

S-P! S-p 2 

We show in another section how any strictly proper rational function (that is, the polynomial 
order in the numerator is less than that of the denominator) in s can be put into this form. 
The connection between Cauchy's theorem and the residue theorem is indicated in the figure 
below. Here, the integral inside the closed path is zero because it excludes the three simple 
poles. Without the cuts however, the integral over the outermost curve would be equal to 
the summation of residues for the poles within (multiplied by i2n). Note that it is only the 
terms with (s — a^) -1 , i.e., simple poles at <2j, that generate nonzero components. 



22.4 Background for the Mapping Theorem 
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Im(s) 




Re(s) 



Looking toward the mapping theorem now, consider the function 



(s- Zl ) m ^{s- z 2 ) m * ... 



(s- Pl ) n i(s-p 2 ) n * ... ' 
Working through some algebra, it is straightforward that 



/'(«)//(*) 



mi 



+ 



m 2 



+ ... - 



Tlx 



n 2 



S — Z\ s — z 2 



s - pi s - p 2 



resulting in 



J f(s)/f(s)ds = i27r(Z-P), 



where Z = m\ + m 2 + . . . is the number of zeros in /, and P = n\ + n 2 + . . . is the number 
of zeros. The mapping theorem comes from now looking in detail at the integral above: 



m 
m 

f'(s)/f(s) 



\f(s)\e^ 

4M.«M+i|/( a )| e ««^ 



ds 

1 d\f(s)\ dd(s) 

~T~ ^ ~ 



ds 



(is ^ 

dlog |/( a )| d»(a) 
+ z- 



ds 



Considering the integral of f'(s)/f(s) over a closed contour in the s-plane, we take advantage 
of the fact that exact differentials rflog \f(s) \ and d9(s) have been found. Both terms pertain 
to the f(s) plane, not the f'(s)/f(s) plane. The first integral is zero, 
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Tdlog |/|=0, 

J Sl 

since log\f(s) \ has to be the same at the start and end of the closed contour. 
Taking 6(s) counterclockwise through angle u2tt results in the second term 

rn2n 

I id6 = in2n. 
Jo 

As noted above, a single circuit in the s-plane may or may not map to a single circuit in the 
/(s)-plane, so n depends directly on the function /(s), and is not necessarily one. Assembling 
the results, we have 

i2n(Z -P) = in2n — > Z - P = n, 

which is the mapping theorem. In words, it states that the number of zeros minus the number 
of poles enclosed in a contour in the s-plane is equal to the number of counterclockwise 
encirclements of the origin, in the f(s) plane. 

23 APPENDIX 2: ADDED MASS VIA LAGRANGIAN 
DYNAMICS 

The development of rigid body inertial dynamics presented in a previous section depends on 
the rates of change of vectors expressed in a moving frame, specifically that of the vehicle. An 
alternative approach is to use the lagrangian, wherein the dynamic behavior follows directly 
from consideration of the kinetic co-energy of the vehicle; the end result is exactly the same. 
Since the body dynamics were already developed, we here develop the lagrangian technique, 
using the analogous example of added mass terms. Among other effects, the equations elicit 
the origins of the Munk moment. 

23.1 Kinetic Energy of the Fluid 

The added mass matrix for a body in six degrees of freedom is expressed as the matrix M a , 
whose negative is equal to: 





Xy 




X p 




Xf 


Yu 


n 


Y- 

1 w 


Yi> 




Yf 




Zy 




Z p 


z* 




Ku 


K t 




K p 






Mu 


My 




M p 




Mr 




N v 




N p 


N 4 


N r 



(260) 



23. 1 Kinetic Energy of the Fluid 
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where (X, Y, Z) denotes the force, (K, M, N) the moment, (u, v, w) denotes the velocity and 
(p, q, r) the angular velocity. The sense of M a is that the fluid forces due to added mass are 
given by 



( Y } 

^-am 




U 


Y 

1 am 




V 




„ d 
• = - M "dt ' 


w 

> 


K 


V 


M am 




q 


N 

K_ 1 v am J 







(261) 



The added mass matrix M a is completely analagous to the actual mass matrix of the vehicle, 





m 











zg 


-Vg ' 







m 





-zq 





x G 


M = 








m 


Vg 


-x G 








-ZG 


Vg 


Ixx 


Ixy 


Ixz 




zg 





-x G 






lyz 




. -Vg 


x G 





Ixz 


lyz 


Izz _ 



(262) 



where [x Gl y Gl z G ] are the (vessel frame) coordinates of the center of gravity. The mass matrix 
is symmetric, nonsingular, and positive definite. These properties are also true for the added 
mass matrix M a , although symmetry fails when there is a constant forward speed. 
The kinetic co-energy of the fluid E k is found as: 



E k = ~q T M a q (263) 
where q T = (u,v,w,p,q,r). We expand to find in the non-symmetric case: 



— 2E k = XuU 2 + Xi,uv + X^uw + Xpup + XgUq + X^ur + 
Yuuv + YyV 2 + Yy,vw + Ypvp + YgVq + Y^vr + 
ZyUW + Z^vw + Z^w 1 + Zpiup + Zqtvq + Z^wr + 
Kuup + Kyvp + K^wp + Kpp 2 + Kqpq + Kj-pr + 
Muuq + Myvq + M^wq + Mppq + M g q 2 + M^qr + 
Nuur + Ni,vr + N^wr + Nppr + N 4 qr + Nf.r 2 . (264) 

For the purposes of this discussion, we will assume from here on that M a is symmetric, that 
is M a = Mj, or, for example, Y^ = Xy. In general, this implies that the z'th force due to 
the fth acceleration is equal to the j'th force due to the i'th acceleration, for % and j = 
[1,2,3,4,5,6]. Refering to the above equation, these terms occur as pairs, so the asymmetric 
case could be reconstructed easily from what follows. By restricting the added mass to be 
symmetric, then, we find: 
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-2E k = XuU 2 + YyV 2 + ZwW 2 + Kpp 2 + M q q 2 + Nj.r 2 + 
2X i] uv + 2X^uw + 2Xpup + 2X q uq + 2Xf.ur + 
2Y lb vw + 2YpVp + 2Y q vq + 2YfVr + 
2Z p wp + 2ZgWq + 2ZrWr + 
2K q pq + 2i^pr + 

2M+qr. (265) 



23.2 Kirchhoff's Relations 

To derive the fluid inertia terms in the body-referenced equations of motion, we use Kirch- 
hoff's relations with the co-energy E k ; see the derivation below, or Milne-Thomson (1960). 
These relations state that if v — [u,v,w] denotes the velocity vector and uj = [p, q, r] the 
angular velocity vector, then the inertia terms expressed in axes affixed to a moving vehicle 
are 



9E k ,^rr\ 

-fx— , (266) 

Q = -£l^|-axf|-* X ^. (267) 

where F = [X, Y, Z] express the force vector, Q = [K, M, N] the moment vector, and x 
denotes the cross product. 




23.3 Fluid Inertia Terms 

Applying Kirchhoff's relations to the expression for the kinetic co-energy with a symmetric 
added mass matrix, we derive the following terms containing the fluid inertia forces: 



X = +XuU + Xi(v - ur) + Xu,(w + uq) + X p p + X q q + Xff 
—YyVr + Yw{vq — wr) — Yppr — Y q qr — Y^r 2 

+Z lb wq + Zppq + Z q q 2 + Z^qr. (268) 
The Y and Z forces can be obtained, through rotational symmetry, in the form: 



Y = +X v (u + ur) - Xyj)up 

+Y v {v + vr) + Y&(w - vp + wr) + Y p (p + pr) + Y q (q + qr) + Y+\ 
-Z^wp- Z p p 2 - Z q pq - Zrpr 



(269) 



23.4 Derivation of Kirchhoff's Relations 
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Z = —X^uq + Xv(up — vq) + X&(u — wq) — Xppq — X q q 2 — X^qr 
+Yi,vp + Yyj{i) + wp) + Ypp 2 + Ygpq + Y+pr 

+ZyjW + Zpp + Zg<j + Zri- (270) 

The apparent imbalance of coefficients comes from symmetry, which allows us to use only 
the 21 upper-right elements of the added mass matrix in Equation 260, e.g., M{, = Yq. 
The moments K, M, and N are obtained in a similar manner as: 



K = —X^wu + X^uv + Xj.uq + Xpii — XqUr 

—Yi,vw + Yw(v 2 — w 2 ) + Yp(v — wp) — Y q (vr + wq) + Yf(vq — wr) 
+Z lb vw + Zp{w + vp) + Z q {vq — wr) + Zf.(wq + vr) 
+K p p + K q (q- rp) + Kf(f + pq) 
+Mr{q 2 - r 2 ) - M q qr 

+Nf.qr (271) 



M = +XuUW + Xi,vw + Xu,(w 2 — u 2 ) + Xp(ur + wp) + X q {u + wq) + Xr(wr — up) 
—Y^uv + Ypvr + Yqi) — Y^vp 

—ZwUW + Z v [wr — up) + Z q {w — uq) — Zf-(ur + wp) 
+Kppr + K q (p + qr) + Kf(r 2 - p 2 ) 
+M q q - Mfpq + M+r 

-Nfpr (272) 



N = —XuUV + Xi,{u 2 — v 2 ) — X^vw — Xp{uq + vp) + X q {up — vq) + Xf{u — vr) 
+Y i) uv + Y^uw + Yp(up - vq) + Y q (uq + vp) + Y+{v + ur) 
—Zpwq + ZgWp + ZfW 
-K m + Kq{p 2 -q 2 ) + Kr(p- qr) 
+Mq P q + Mr(q + pr) 

+N+f (273) 



23.4 Derivation of Kirchhoff's Relations 

We can derive Kirchhoff's relation for a lagrangian L(v,u,t), involving the velocity v and 
angular velocity uj, whose components will be expressed in a local coordinate system rotating 
with the angular velocity uj, i.e. in a reference system fixed on the body. The principle to 
satisfy is that of least action, i.e. to minimize the integral / (Crandall et a/., 1968): 
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1= [** L(v,u,t)dt (274) 

At the minimum value of / - the admissible condition - the variation of /, SI, and hence of 
L with the velocity v and angular rate u vanishes. 

Our condition SI = 0, as written, involves only the lagrangian, which in the more general case 
is the kinetic co-energy minus the potential energy of the system. Since we are considering 
the motion of a body in an unbounded, homogeneous fluid, there is no potential energy, so 
the lagrangian is exactly the kinetic energy: 

L = E k . (275) 

Hamilton's principle in its general form also accounts for applied forces and moments S. They 
are defined to align with the generalized, infinitesimal linear and angular displacements Sfj 
and S(f), leading to 



SI — f SL(v, lj, t) + (H U)t))W , Sfj) + S(p 



dt. 



(276) 



Now, the lagrangian is invariant under coordinate transformation, so it is a function of the 
free vectors of velocity and angular velocity. Using the notation detailed in the Nomen- 
clature section below, Sfj and S<j) are interpreted as free vectors, while Srj and 5(f> are the 
projections of the free vectors onto a given reference frame. The following relationships link 
the displacements with the body-referenced rates: 



v = 



dt] 

■7= + u <g) rj, and 
at — 

~dt' 



(277) 
(278) 



A variation of the lagrangian at a given time t is, to first order, 



SL(v^t) = ^Sv\ + /^Su 



(279) 



The variations Sv and Su, in view of equations 277 and 278, can be written as 



Sv 



StFx + {iFSx} 



' dSrj 
~dt 



+ u <S> Srj I x + v x S(f>, 



(280) 
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Su = 5uFx + {uF5x} 

(d56\ T 
= (^=j x + ux56. 



(281) 



The terms {rFSx} and {uFSx} above represent the effects of the variation of body orienta- 
tion, and lead to one of the more subtle points of the derivation. It can be shown that 5x, 
the displacement of the unit triad x is 56 x x, leading for example to v T (56 x x). This form, 
however, fails to capture what is meant by {iF5x} and {uFSx}, specifically: how the free 
vectors 5v and 5w change as the triad rotates, but the projections y_ and w_ remain constant. 
With this in mind, one can easily derive the proper interpretations, as written above. 
We now return to the evaluation of the lagrangian L. Combining terms, we see that there 
will be five inner products to consider. Writing the terms involving the time derivatives as 
5Li and 5Ii, we have from an integration by parts 



Sh 




8L fdSrfy \ + /8L fdSf 



dv\dt J 



diJ \ dt 



x 



dt 



ddL \ I ddL _ T \ 



dt. 



(282) 



There are no terms remaining at the time boundaries because the lagrangian is zero at these 
points. 

In evaluating the remaining terms, in 5I 2 , one needs to use the following property of triple 
vector products: 



< a,b x c >=< c,b x a >=< b,cx a > . 
Hence, the variation of the action I 2 , expressed in terms of vectors, is 



(283) 



51, = 



*2 



tl 



rt2 
Jti 



x 5rj + v x 5^ + x 5^ dt 

v x — ,5r] j + I v x — + u x — ,56) dt. 



(284) 



Finally, we write the part of 51 due to the generalized forces: 




(285) 



The Kirchoff relations follow directly from combining the three action variations. The kinetic 
co-energy we developed is that of the fluid, and so the generalized forces indicate what forces 
the body would exert on the fluid; with a sign change, these are the forces that the fluid 
exerts on the body. 
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23.5 Nomenclature 

23.5.1 Free versus Column Vector 

We make the distinction between a free vector /, which is an element of a linear vector 
field, and a column vector / which denotes the components of the vector / in a given 
coordinate system. The connection between the two concepts is given via the free triad x, 
containing as elements the unit vectors of the chosen Cartesian system, k, i.e.: 



x = (i,j,k) T (286) 
The notation is hybrid, but convenient since it allows us to write: 



/ = fx (287) 



where / T denotes the transpose of /, i.e. a row vector. The product between the row vector 
f T and the column vector x is in the usual matrix multiplication sense. 

23.5.2 Derivative of a Scalar with Respect to a Vector 

The derivative of a scalar L with respect to a vector x is denoted as: 

dL , s 

m (288) 

and is a vector with the same dimension as x, whose element i is the derivative of L with 
respect to the ith element of x, and in the same direction: 

fdL\ dL_ 
\dx J . dxi 

23.5.3 Dot and Cross Product 

— * — * 

The dot product of two free vectors / and g is denoted as < /, g >, and we use the following 
notation for column vectors: < / T f,£ T x >. The cross product of two vectors / and g is 

z ' = f x 9i denoted in terms of column vectors as z = f <S> g- 

24 APPENDIX 3: LQR VIA DYNAMIC PROGRAM- 
MING 



There are at least two conventional derivations for the LQR; we present here one based on 
dynamic programming, due to R. Bellman. The key observation is best given through a 
loose example. 



24. 1 Example in the Case of Discrete States 
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24.1 Example in the Case of Discrete States 

Suppose that we are driving from Point A to Point C, and we ask what is the shortest path 
in miles. If A and C represent Los Angeles and Boston, for example, there are many paths 
to choose from! Assume that one way or another we have found the best path, and that a 
Point B lies along this path, say Las Vegas. Let X be an arbitrary point east of Las Vegas. 
If we were to now solve the optimization problem for getting from only Las Vegas to Boston, 
this same arbitrary point X would be along the new optimal path as well. 
The point is a subtle one: the optimization problem from Las Vegas to Boston is easier than 
that from Los Angeles to Boston, and the idea is to use this property backwards through 
time to evolve the optimal path, beginning in Boston. 

Example: Nodal Travel. We now add some structure to the above experiment. Consider 
now traveling from point A (Los Angeles) to Point D (Boston). Suppose there are only three 
places to cross the Rocky Mountains, B\, B 2 , B%, and three places to cross the Mississippi 
River, Ci, C 2 , C3. 3 By way of notation, we say that the path from A to B\ is AB\. Suppose 
that all of the paths (and distances) from A to the 5-nodes are known, as are those from the 
.B-nodes to the C-nodes, and the C-nodes to the terminal point D. There are nine unique 
paths from A to D. 

A brute-force approach sums up the total distance for all the possible paths, and picks the 
shortest one. In terms of computations, we could summarize that this method requires nine 
additions of three numbers, equivalent to eighteen additions of two numbers. The comparison 
of numbers is relatively cheap. 

The dynamic programming approach has two steps. First, from each 5- node, pick the best 
path to D. There are three possible paths from Bi to D, for example, and nine paths total 
from the B-\eve\ to D. Store the best paths as BiD\ opt , B 2 D\ opt , B 3 D\ opt . This operation 
involves nine additions of two numbers. 

Second, compute the distance for each of the possible paths from A to D, constrained to the 
optimal paths from the B-nodes onward: ABi + BiD\ opt , AB 2 + B 2 D\ opt , or AB 3 + B 3 D\ opt . 
The combined path with the shortest distance is the total solution; this second step involves 
three sums of two numbers, and total optimization is done in twelve additions of two numbers. 
Needless to say, this example gives only a mild advantage to the dynamic programming 



3 Apologies to readers not familiar with American geography. 
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approach over brute force. The gap widens vastly, however, as one increases the dimensions of 
the solution space. In general, if there are s layers of nodes (e.g., rivers or mountain ranges), 
and each has width n (e.g., n river crossing points), the brute force approach will take (sn s ) 
additions, while the dynamic programming procedure involves only (n 2 (s — 1) + n) additions. 
In the case of n = 5, s = 5, brute force requires 6250 additions; dynamic programming needs 
only 105. 

24.2 Dynamic Programming and Full-State Feedback 

We consider here the regulation problem, that is, of keeping x desired = 0. The closed-loop 
system thus is intended to reject disturbances and recover from initial conditions, but not 
necessarily follow y-trajectories. There are several necessary definitions. First we define an 
instantaneous penalty function l(x(t),u(t)), which is to be greater than zero for all nonzero 
x and u. The cost associated with this penalty, along an optimal trajectory, is 

POO 

J= / l(x(t),u(t))dt, (290) 
Jo 

i.e., the integral over time of the instantaneous penalty. Finally, the optimal return is the 
cost of the optimal trajectory remaining after time t: 

POO 

V(x(t),u(t)) = J 1(x(t),u(t))cIt. (291) 

We have directly from the dynamic programming principle 

V(x(t),u(t)) = mm {l(x(t),u(t))St + V(x(t + 5t),u(t + St))} . (292) 

The minimization of V(x(t),u(t)) is made by considering all the possible control inputs u in 
the time interval (t, t + St). As suggested by dynamic programming, the return at time t is 
constructed from the return at t + St, and the differential component due to l(x,u). If V is 
smooth and has no explicit dependence on t, as written, then 

f)V (It 

V(x(t + St),u(t + St)) = V(x(t),u(t)) + — — St + h.o.t. — > (293) 

\J JL (JjIj 

dV 

= V(x(t),u(t)) + -z-(Ax(t) + Bu(t))St. 
ox 

Now control input u in the interval (t,t + St) cannot affect V(x(t),u(t)), so inserting the 
above and making a cancellation gives 

= mm ^l(x(t),u(t)) + ^(Ax(t) + Bu(t))^ . (294) 
We next make the assumption that V(x,u) has the following form: 

V(x,u) = \x T Px, (295) 
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where P is a symmetric matrix, and positive definite. 45 It follows that 

% = *P^ (296) 



= mm {l(x,u) + x T P (Ax + Bu)} 



We finally specify the instantaneous penalty function. The LQR employs the special quadratic 
form 

l(x, u) = -x T Qx + )-u T Ru, (297) 

where Q and R are both symmetric and positive definite. The matrices Q and R are to be 
set by the user, and represent the main "tuning knobs" for the LQR. Substitution of this 
form into the above equation, and setting the derivative with respect to u to zero gives 



= u T R + x T PB (298) 
u T = —x T PBR~ 1 
u = -R- 1 B T Px. 

The gain matrix for the feedback control is thus K = R~ l B T P. Inserting this solution 
back into equation 297, and eliminating u in favor of x, we have 

= -x T Qx - -x T PBR- 1 B T P + x T PAx. 
2 ^ 2 

All the matrices here are symmetric except for PA; since x T PAx = x T A T Px, we can make 
its effect symmetric by letting 

x T PAx = -x T PAx + -x T A T Px, 
2 2 

leading to the final matrix-only result 

= Q + PA + A T P - PBR- l B T P. (299) 



25 Further Robustness of the LQR 

The most common robustness measures attributed to the LQR are a one-half gain reduction 
in any input channel, an infinite gain amplification in any input channel, or a phase error 
of plus or minus sixty degrees in any input channel. While these are general properties 

4 Positive definiteness means that x T Px > for all nonzero x, and x T Px = if x = 0. 

5 This suggested form for the optimal return can be confirmed after the optimal controller is derived. 
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that have a clear graphical implication on the Bode or Nyquist plot, other useful robustness 
conditions can be developed. These include robustness to uncertainty in the real coefficients 
of the model (e.g., coefficients in the A matrix), and certain nonlinearities, including control 
switching and saturation. We will use the Lyapunov stability and the LM1 formulation of 
matrix problems in this section to expand these ideas. 

Saturation nonlinearities in particular are ubiquitous in control system application; we find 
them in "railed" conditions of amplifiers, rate and position limits in control surface actuators, 
and in well-meaning but ad hoc software limits. As shown below, moderate robustness in 
saturation follows from the basic analysis, but much stronger results can be obtained with 
new tools. 

When the LQR is used to define the loop shape in the loop transfer recovery method (as 
opposed to the Kalman filter in the more common case), these robustness measures hold. 

25.1 Tools 

25.1.1 Lyapunov's Second Method 

The idea of Lyapunov's Second Method is the following: if a positive definite function of the 
state x can be found, with V(x) = only when x — 0, and if dV{x)/dt < for all time, then 
the system is stable. A useful analogy for the Lyapunov function V(x) is energy, dissipated in 
a passive mechanical system by damping, or in a passive electrical system through resistance. 

25.1.2 Matrix Inequality Definition 

Inequalities in the case of matrices are in the sense of positive and negative (semi) definite- 
ness. Positive definite A means x 7 Ax > for all x; positive semidefinite A means x T Ax > 
for all x. With A and B square and of the same dimension, 

A < B means x T Ax < x T Bx, for all x. (300) 
Also, we say for the case of a scalar 7, 

A < 7 means A - 7/ < 0. (301) 

25.1.3 Franklin Inequality 

A theorem we can use to assist in the Lyapunov analysis is the following, attributed to 
Franklin (1969). 

A T B + B T A < jA T A + -B T B, for all real 7 > 0. (302) 

7 

The scalar 7 is a free parameter we can specify. It is assumed that the matrices A and B 
are of consistent dimensions. 



25.1 Tools 
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25.1.4 Schur Complement 

Consider the symmetric block matrix defined as follows: 



M = 



A B 
B T D 



(303) 



The Schur complements of blocks A and D are defined, respectively, as 



D - B T A- 1 B 
A — BD^B T , 



and an important property is that 

M > i — > T A > and r n > and A > and D > 0. 



(304) 



(305) 



The Schur complements thus capture the sign of M, but in a smaller dimensioned matrix 
equation. The fact that both A and D have to be positive definite, independent of the Schur 
complements, is obvious by considering x that involve only A or D, e.g., x = [lixn A 0i x „ D ] T , 
where ua and rip are the dimensions of A and D. 

25.1.5 Proof of Schur Complement Sign 

It is easy to verify that 



M 



A 


B ' 




B T 


D 





I 

B T A~ l I 



A 

D-B T A- 1 B 



I A- l B 
/ 



The outer matrices are nonsingular so there is a one-to-one transformation 



/ A- l B 
I 



x, 



(306) 



(307) 



and thus M > is equivalent to A > and T A > 0. The proof for T D is completely 
analogous. 

25.1.6 Schur Complement of a Nine-Block Matrix 

For the purposes of derivation below, consider now the symmetric matrix 



M 



ABC 
B T D 
C T E 



(308) 



Using the Schur complement of the diagonal block including D and E, we have 

M > < — > A - BD- l B T - CE^C 7 > and D > and E > 0. (309) 
Again, positive definiteness of a large matrix is reflected in a smaller matrix inequality. 



136 



25 FURTHER ROBUSTNESS OF THE LQR 



25.1.7 Quadratic Optimization with a Linear Constraint 

Consider the quadratic function V = x T Px, P symmetric and positive definite, and the 
minimization of V subject to a constraint c T x = a, where c T is a row vector, and a is a 
scalar. The solution is: 

min V = (310) 
To show that this is true, augment the objective with a LaGrange multiplier to become 

V = V + \{c T x-a). (311) 

Setting dV/dx = gives 2x T P + Xc T = 0. Post-multiply this equation by P _1 c, and solve 
to give A = —2a/c T P~ 1 c. Then, using 2x T P + Xc T = again, post-multiply by x/2, and 
insert A, concluding the proof. 



25.2 Comments on Linear Matrix Inequalities (LMI's) 

The reason that solving M > is simpler than the equivalent Schur complement inequality 
is that, as long as M{y) is affine in the unknowns y, the solution set is convex. Affine means 
that 



M(y) = M + K[yy...y\N, (312) 

i.e., that M is a constant matrix plus terms linear in each element of y. Note the matrices 
K and N here, as well as the columnwise expansion of y, are needed to put elements of y in 
arbitrary locations of M. 

Convex means that if y x and y 2 are feasible solutions to the inequality, then so is y 3 = 
ayi + (1 — Oi)y 2) for any a between zero and one. A convex solution set has no hidden 
corners or shadows, because every point on a straight line between two solutions is also 
a solution. Computer programs solving convex optimization problems can always get the 
global optimum efficiently, or determine that no solution exists. 

From the Schur complements above, clearly we can transform a (suitable) matrix equation 
that is quadratic or higher in the unknowns into an affine M(y), and therefore make an 
efficient solution. The term linear matrix inequality refers of course to the fact that M is 
affine and hence linear in the unknown variables, but also to the methods of analysis and 
numerical solution, wherein the idea is to recognize 



M(y) = M + M m + M 2 y 2 + ...M n y n , 



(313) 



where n is the dimension of the unknown vector y. 
Example. Consider the LMI 



M(y) 



Vi 

2/2 



1 27/ 2 -l 

i V2-yi 



> o. 



(314) 
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We use the theorem that M > is achieved if each of the submatrices cornered at the (1,1) 
element has positive determinant. Then M > is equivalent to 



yi > i 

- 1)0/2 -?/i)-(2y 2 - l)0/ 2 + l) > 0. (315) 

The solution set is shown as S in the figure below; M is affine in y and the solution set is 
convex. 




25.3 Parametric Uncertainty in A and B Matrices 
25.3.1 General Case 

Let S be the solution to the matrix Riccati equation 

A T S + SA + Q- PBR- l B T P = 0. (316) 

We consider here perturbations to the design plant matrices A and B of the form x = 
(A + AA)x + (B + AB)u. The control in the LQR is u = -R- l B T Px. Select a Lyapunov 
function of the form V(x) = x T Sx. Then 



V — x Sx 



->T 

X 



A T S + SA + AA T S + SAA - 2SBR- 1 B T S- 



SBR- 1 AB T S - SABR^BS 



x. 



(317) 



We give the following special structures to the perturbations: 
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AA = D A A A E A 
AB = D B A B E B . 



(318) 



The idea is to capture all of the uncertainty in and A B and capture the structure of the 
perturbations using the other four matrices. We have, using the Franklin inequality, 



V = x 



< x T 



A'S + SA- 2SBR~ 1 B I S + E A A A D A S + SD A A A E A 
-SBR^ElAlDlS - SDbAbEbR^BS 



(319) 



x 



1 



A T S + SA- 2SBR~ 1 B T S + lA E T A A T A A A E A + —SD A D T A S 

lA 



+~i B SBR- 1 E T B A T B A B E B R- l B T S + —SD B D T B S 

IB 



X 



< x T 



1 



A 1 S + SA - 2SBR- L B 1 S + j A E A E A + — SD A D A S 

lA 



+~1 B SBR- 1 E T B E B R-'B 1 S + —SD B D B S 

IB 



-1 dT, 



1 



X, 



(320) 



where the last inequality holds if 



A T A A A < L 



A A ^ 



ArAr < /. 



(321) 
(322) 



For a given perturbation model with matrices D A , E A , D B , and E B , we thus obtain a matrix 
inequality which, if met, guarantees stability of the system for all x: 



A T S + SA- 2SBR- 1 B T S + y A E A E A + —SD A D T A S 

lA 

+ r y B SBR~ 1 E B E B R~ 1 B T S + —SD B D T B S < 0. 

IB 



(323) 



Since scalars ^ A and j B can take on any positive value, they remain to be picked, and in 
fact judicious choices are needed to maximize the robustness guarantees - i.e., to minimize 
the conservativeness. 



25.3.2 Uncertainty in B 

Putting aside A A for the moment, and substituting in the Ricatti equation into the above 
leads to 



25.3 Parametric Uncertainty in A and B Matrices 
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-Q - SBR- l B T S + -i B SBR- l ElE B R- 1 B T S + —SD B D T B S < 0. (324) 

Following the presentation of gain and phase margins of the LQR in a previous section, we 
here consider gain margins through the diagonal matrix N, such that 

A B = BN = By/R-^Nl x I x yj R\N\ = D B A B E B . (325) 

\N\ denotes the matrix made up of the absolute values of the elements of N. This A B is 
a very specific structure, whose rationale becomes evident below. The stability inequality 
becomes 

-Q - SBR- l B T S + - fB SB\N\R- 1 B T S + —SBR'^N^S < 0, (326) 

IB 

or, equivalently, 

/ - 1b\N\ |JV| > 0, (327) 

7b 

because Q can be arbitrarily small in the LQR; certainly a large value of Q increases the 
robustness (see the case of A A below, for example). Since N is diagonal, it is sufficient to 
keep all 

1-|AT M | (lB + ^j > 0. (328) 

Solving for j B gives 



l±Jl- 4iV2. 

" B = 2 A, • ^ 329 ) 

showing that real, positive j B are attained for | JV^^ | < 1/2. Hence, this analysis confirms 
the one-half gain reduction property derived earlier. On the other hand, it confers only a 
one-half gain amplification robustness, which is overly conservative. Usage of the Franklin 
inequality to symmetrize the components involving AB is the cause of this, since it can be 
verified that using Equation 319 directly gives 



-Q - SBR-\I + iV + N)B T S < 0, or 

l + 2iV M > 0, (330) 

showing both the one-half reduction and infinite upward gain margins. 

In summary, in the special case of diagonal N in the model of AB, we crafted D B and 
E B so as to align terms of the form SBR~ 1 B T S, simplifying the analysis. This leads to a 
conservative result using the Franklin inequality, but the full, standard result on a second 
look. Better results could be obtained for specific cases, but probably with greater effort. 
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25.3.3 Uncertainty in A 

Uncertainty in the system matrix A has a different role in robustness than does uncertainty 
in B, which is primarily an issue of gain and does not affect the open-loop stability. Pertur- 
bations to the design matrix A cause the open- and closed-loop poles to move, potentially 
destabilizing either or both of the open- and closed-loop systems. Consider the main in- 
equality for stability again, setting aside AB now: 

-Q - SBR- 1 B T S + jaE^Ea + —SD A D T A S < 0. (331) 

lA 

Aligning the uncertain terms with SBR~ l B T S (as with AB above) proves to be impractical 
here because it imposes a specific structure on AA, and requires that elements of AA would 
have to change together. A more useful result involves a diagonal A A , and for the purposes 
of illustration, we consider a second order system as follows: 



D A A A E A = 

di5i d 2 S 2 






d 2 ' 




" 5 1 


" 













s 2 



(332) 
(333) 



In addition, we set B = [b 0] T , so that R is a scalar r, and we set Q diagonal. The resultant 
inequality is 



-Q + j A I-S 



' b 2 /r " 










{d\ + dl)/ lA 




S < 0. 



(334) 



The condition for stability will be met if Q — ^ A I > and b 2 /r — (d\ + d^/jA > 0, or 
equivalently, setting ■ja = m i n Qi,i, if 



b 2 

— min Qi i > d\ + d\. (335) 



Now let, for example, 



B 

Q 
R 



-0.5 - 
1 

[1 Of 
diag(l, 1) 
1. 



-0.5 




This informs us that the robustness condition is d 2 + d 2 < 2, and very substantial perturba- 
tions in A are allowed, for example 
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A + AA 



0.5 
1 



0.5 




Computed results that confirm the bound are shown in Figure 9. An interesting note is that 
the gain reduction margin in the presence of this AA suffers dramatically, with failure of 
the robustness if b < 0.93. Intuitively, the uncertainty in A has used up almost all of the 
available safety margin. 




Figure 9: (left) open-loop pole locations along the boundary of the stability region in Equa- 
tion 335; (right) closed-loop pole locations with LQR. 



25.3.4 A and B Perturbations as an LMI 

Using the triple LMI of Equation 309, we have, for the combined A and B structured 
uncertainty models, 



1 



-Q - SBR-'B 1 S + j A E A E A + —SD A D\S 

lA 

+- fB SBR- 1 ElE B R- 1 B T S + —SD B D T B S < 

IB 



is equivalent to 



(336) 



-Q- SBR^B 1 S + ^ A E 1 A E A + lB SBR~ 1 E B E B R- 1 B 1 S SD A SD B 



D T A S 
DlS 



-IaI 






-J B I 



< 0. (337) 



The solution set [7^ j B ], if it exists, is convex. If a solution does not exist, one could 
lower the robustness by altering D A , E A , D B , or E B . Alternatively, increasing Q, the state 
penalty, where possible may increase robustness. 
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25 FURTHER ROBUSTNESS OF THE LQR 



25.4 Input Nonlinearities 

Control systems in practice often make use of switched or saturated action in the control 
channels. Switched and saturated control channels are defined respectively as 



switch(n, U) = (7sign(tt), 



(338) 



and 



sat(u, U) = Us&t(u/U, 1) = I • / \ " 
v ' v ' ' [ (/sign (it), \u 



u\<U 
> U. ' 



(339) 



where the scalar U is positive. If the second argument of sat(-) is omitted, it is taken to be 
one. 



u 


switch(u,U) 






u 

-U 





U 

-u 


i sat(u,U) 










u 
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Our stability analysis so far is one wherein we seek to guarantee a negative rate in V(x) for 
all time. Referencing Equation 330, it is clear that any N can be tolerated within the range 
(—.5,oo), even if it is time varying. It follows directly then that switching and saturated 
control channels cannot destabilize the LQR system until the gain reduction is more than 
one-half, that is, \u\ > 2(7. In practice, however, systems routinely operate beyond this point, 
and this is the aspect of saturated control in particular that we wish to pursue now. This 
section details the approach of Huang & Lam (2002), which uses LMI's and the structured 
uncertainty models above. 

We say that f(u) is the diagonal function taking control commands from the LQR (u) to 
the input to the plant; elements of / correspond to the saturated operator above. Through 
scaling the B matrix, we can set Ui = 1, and hence /'s extreme outputs are ±1. 
First, we assert that under the proper choice of positive definite matrix W, and for certain 
vectors z, the following inequality holds: 



2z T f{R' 1 z) > z T WR' l z. 



(340) 



If the control penalty matrix R is diagonal, and we constrain W also to be diagonal, this 
condition is met if, for each control channel, 



2z ? 'Sat 



Zi 



R 



> z 



. Wi,i 

i id ' 



(341) 
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If the channel is not saturated, then we require W i:i < 2. If the channel is saturated, then it 
is required instead that 

> Wi,i. (342) 

We will use the inequality of Equation 340 in our Lyapunov analysis, with special attention 
to these conditions. Because for the LQR we have Zi = e[B T Sx, where is the unit vector 
in the i'th direction, Equation 342 effectively places a bound on the state under saturated 
control, namely 

\ejB T Sx\ < 2R i jW i , i . (343) 

Viewed this way, clearly we would like to make each W i:i as small as possible. As we show 
below, however, the matrix W appears in the only negative term of V(x); the optimization 
then is to maximize the state space satisfying Equation 342, while keeping W big enough to 
guarantee stability. 

With V(x) = x T Sx, and considering uncertainty AB only, we have 

V{x) = x T {SA + A T S)x-2x T SBf{R- 1 B T Sx)-2x T SABf{R- 1 B T Sx). (344) 

Note that AB never appears inside / because / operates on the control it, calculated with 
B alone. Applying the Franklin inequality to the last term yields 

V{x) < x T {SA + A T S)x-2x T SBf{R- 1 B T Sx) 

+j B x T SD B DlSx + —f{R~ l B T Sx)ElE B f{R' l B T Sx)). (345) 

7b 

Next, consider the remaining asymmetric term. Applying Equation 340, under the proper 
constraints on W, we obtain 

V(x) < x T {SA + A T S)x-x T SBWR' 1 B T Sx 

+j B x T SD B DlSx + —f{R- l B T Sx)ElE B f(R- l B T Sx)). (346) 

7b 

Finally, the last term is bounded, giving 

SA + A T S - SBWR- 1 B T S + 7 B SD B D T B S + —SBR- l E T B E B R- l B T S < (347) 

7b 

as the sufficient condition for stability. This last step is extremely conservative, since in fact 
\fi\ < 1; Huang and Lam address this point. 
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The supporting conditions on elements of W are linear: W i;i must be as in Equation 342, 
and less than two. As noted above, these translate to linear constraints on x. Call the set 
of states that satisfy these conditions fy. The problem now is that although the Lyapuonov 
function derivative is clearly less than zero for states within the state trajectory could 
leave * . This is entirely likely since V is a quadratic function, leading to elliptical contours 
as in the figure below. 



V contours 




As a result, the range of states leading to guaranteed stability is limited to those that are 
known not to leave In other words, we have to settle for the largest ellipse completely 
contained in ty. We now invoke Equation 310, referencing Equation 343 as the set of linear 
constraints, to yield 



AR Z - 

V ™ = ^ W^BTSBe, (348) 

We would like to maximize V max , which is equivalent to maximizing the minimum value 
of all the possible right hand sides of the above equation. This itself can be posed as an 
optimization in a new scalar k: 

minimize k such that kl - )-R~ l WT > 0, (349) 

where T is diagonal, with T i:i = \j ej B T SBei. It follows that any x is stable when V(x) = 
x T Sx<l/k 2 . 

In summary, we have three constraints and one optimization: Equation 347 for Lyapunov 
stability, < W i:i < 2 for cases where channels are not saturated, and Equation 349 for 
maximizing the level of V - an elliptical region in the state space - for which stability is 
guaranteed. 

As in the case of no input saturation, this can be posed and solved as set of LMI's (which 
itself can be posed as a single LMI): 

SA + A T S - SBWR' 1 B T S + lB SD B D T B S SBR~ l El 
E B R- 1 B T S ~i B I 

-R~ X W 
2 



< 


o, 


(350) 


< 


k, 


(351) 


< 


2, 


(352) 


< 


W. 


(353) 
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To solve, pick a value for k and find a feasible solution. Decrease k as far as possible; the 
feasible solution achieves the maximum V from which the range of state space can be found 
fmmV = l/k 2 . 



